Appell-Gibbs Approach in Dynamics of Non-Holonomic Systems

Abstract

Hamiltonian functional and relevant Lagrange’s equations are popular tools in the investigation of dynamic systems. Various generalizations enable to extend the class of problems concerned slightly beyond conventional limits of Hamiltonian system. This strategy is very effective, particularly concerning two-dimensional (2D) and simpler three-dimensional (3D) systems. However, the governing differential systems of most non-holonomic 3D systems suffer from inadequate complexity, when deduced using this way. Any analytical investigation of such a governing system is rather impossible and its physical interpretation can be multivalent. For easier analysis, particularly of systems with non-holonomic constraints, the Appell-Gibbs approach seems to be more effective providing more transparent governing systems. In general, the Appell-Gibbs approach follows from the Gaussian fifth form of the basic principle of dynamics. In this chapter, both Lagrangian and Appell-Gibbs procedures are shortly characterized and later their effectiveness compared on a particular dynamic system of a ball moving inside a spherical cavity under external excitation. Strengths and shortcomings of both procedures are evaluated with respect to applications.

1. Introduction

The energy contained in a dynamic system is given by a scalar potential $\mathcal{E}\left(t\right)$. It is a function of time and system response components (displacement, velocity, and acceleration vectors). Moreover, $\mathcal{E}\left(t\right)$ is a function of system parameters, position in a field of forces (potential or not), internal sources of energy and of the system evolution including a residual energy. The total energy of the system increases or decreases accordingly with external excitation and dissipation of energy. The form of energy contained within the system can have a deterministic or stochastic character and similarly also excitation and dissipation.

Considering the mechanical energy only, the total energy increase/decrease of the system with respect to time should be in equilibrium with the energy supplies and energy losses due to dissipation. This relation can be outlined by the following equilibrium:

$$\frac{\mathrm{d}}{\mathrm{d}t}\left\{\mathcal{E}\left(t\right)\right\}=\mathcal{P}\left(t\right)+\mathcal{S}\left(t\right),$$E1

where $\mathcal{P}\left(t\right)$ is power supply (excitation energy per unity time) and $\mathcal{S}\left(t\right)$ the specific dissipation of energy also per unity time (supposed to be independent on accelerations $\ddot{\mathbf{x}}$). Functions $\mathcal{P}\left(t\right),\mathcal{S}\left(t\right)$ can dispose in special cases with a superior potential, which, however, cannot be incorporated into the potential part of total energy. Eq. (1) has a scalar character.

The energy is a primary value characterizing the system state and its evolution in time. The function $\mathcal{E}\left(t\right)$ and external influences are a background for the derivation of a governing differential system characterizing the system response with respect to initial and boundary conditions. The governing differential system is then deduced from the equivalence of Eq. (1) type using an adequate variational principle. It claims that the form of the system response corresponds with the minimum of energy spent among all admissible shapes of the system reaction. Take a note that many important settings of external forces and dissipation mechanisms do not admit the formulation by means of potentials. In such cases, they should be incorporated separately into the governing differential system using complementary principles and theorems, for example, virtual works, and so on.

We can find in monographs, for example, [1, 2, 3, 4, 5] and many others, various formulations of potentials $\mathcal{E}\left(t\right)$ and functions $\mathcal{P}\left(t\right),\mathcal{S}\left(t\right)$ combining the system parameters (physical and geometric) and the system response vectors $\mathbf{x}$-displacements, $\dot{\mathbf{x}}$-velocities, and $\ddot{\mathbf{x}}$-accelerations. They can be selected in individual cases with respect to physical or geometric complexity of the system, components of the response, which are to be found, deterministic or stochastic character of the system and its excitation, and so on.

2. Basic considerations

Approaches commonly applied to construct mathematical models of dynamic systems with multiple degrees of freedoms (MDOF) follow mostly from principles symbolically outlined by Eq. (1). The equation of this type can be deduced using, for instance, a procedure of virtual displacements. They balance the energy flow in every step and subsequently applied minimization steps try to select such response trajectories, which represent a minimum of energy consumption among all admissible shapes. Let us get briefly through Lagrangian and Appell-Gibbs procedures in order to compare their basic properties. Later, we recognize that most of these properties can be regarded as positive or negative in dependence on a particular problem. Therefore, the solution method should be selected in every particular case very sensitively.

Let us remember that the aim of this chapter is a comparison of Lagrangian and Appell-Gibbs approaches effectiveness to process dynamic systems in holonomic and non-holonomic settings and to help estimate which one is more suitable to be employed in a particular case. Despite that the most important features of non-holonomic systems themselves are briefly treated as well, but for thorough evaluation of their properties, special literature should be addressed. Except five monographs cited in introductory section containing a large number of additional relevant references, a vast number of papers have been published concerning the investigation of various properties of non-holonomic systems.

Motion of an MDOF system with $n$ degrees of freedom can be described by a system of $n$ differential equations and $l$ constraints:

$$\begin{array}{lllll}& m_s{\ddot{x}}_s=X_s+\sum _{r=1}^l{\lambda }_r{A}_{\mathit{rs}},& s=1,..,n,& \mathbf{x}=\left[x_s\right],\mathbf{X}=\mid X_s\mid ,\mathbf{x},\mathbf{X}\in {\mathbb{R}}^n,& \left(a\right)\\ & \sum _{s=1}^n{A}_{\mathit{rs}}{\dot{x}}_s+B_r=0,& r=1,..,l,& \mathit{\lambda }=\left[{\lambda }_r\right],\mathbf{B}=\left[B_r\right],\mathit{\lambda },\mathbf{B}\in {\mathbb{R}}^l,& \left(b\right)\\ & & & \mathbf{A}=\left[A_{\mathit{rs}}\right],\mathbf{A}\in {\mathbb{R}}^{l\times n}.& \end{array}$$E2

Vector $\mathbf{X}$ represents external forces, while $\mathit{\lambda }$ are unknown multipliers. The summation in Eq. (2a) characterizes influence of constraints (holonomic and non-holonomic) related with constraints (Eq. (2b)). These constraints reduce the number of the original degrees of freedom from $n$ to $k=n-l$. The system (Eq. (2)) includes $n+l$ differential equations for $\mathbf{x}$ and $\mathit{\lambda }$ unknown functions $t$, which can be determined, provided $\mathbf{x},\dot{\mathbf{x}}$ are given in an initial point $t_0$. If the system (Eq. (2b)) is fully integrable, it provides $l$ functions $f_r=f_r\left(\mathbf{x},t\right),r=1,..,l$ and constraints can be formulated as $f_r=f_r\left(\mathbf{x},t\right)=c_r$. They are exclusively of a geometric character and the system is holonomic. Corresponding constraints are formulated in displacements only. In principle, $l$ components of $\mathbf{x}$ can be eliminated and then remains to analyze the system with $n-l$ unknowns. Then, it can be considered $\mathit{\lambda }\equiv 0$, and the second part on the right side of Eq. (2a) vanishes. The system with holonomic constraints takes the form:

$$m_s{\ddot{x}}_s=X_s,f_r=f_r\left(\mathbf{x},t\right)=c_r,s=1,..,k,r=1,..,l,k=n-l.$$E3

However, frankly speaking, such an operation is possible rather exceptionally. In general, the full form of Eq. (2) should be treated, despite the system is holonomic. If some (or all) of constraints (Eq. (2b)) are not integrable, then the system is non-holonomic. In practice, we encounter these cases when the formulation of constraints includes velocities (more often velocities only).

We should remember that the non-holonomic constraints introduced in Eq. (2b) represent the most simple version of such constraints, as they are linear and in velocity. Many applications, for example, robotics, wind engineering, automotive systems, plasma physics, and so on, present more complicated types of non-holonomic constraints. Notifications to nonlinear non-holonomic constraints in velocity are given in elderly monographs [2, 3]. Later, many papers have appeared presenting results of systematic research at this field originating from particular physical or engineering problems, for example, [6, 7, 8], where higher derivatives of velocities in non-holonomic constraints are discussed. These attributes have been reflecting also in pure mathematical studies with respect to control theory and systems with a delayed feedback, see, for example, series [9, 10, 11] dealing with generalized Lagrange-d’Alembert-Poincaré equations or other studies devoted to non-holonomic reduction and related problems, see, for example, [12] and many others.

Let us realize now that the virtual work of every constraint force should vanish in the meaning as follows:

$${\lambda }_r\sum _{s=1}^n{A}_{\mathit{rs}}\delta x_s=0,r=1,..,l.$$E4

Therefore, we have with respect to Eq. (2a):

$$\sum _{s=1}^n\left(m_s{\ddot{x}}_s-X_r\right)\delta x_s=0.$$E5

This equation holds for any arbitrary virtual displacements and represents a generalization of the principle of virtual works in statics and of the d’Alembert principle. The important issue is that it does not include any reactions of constraints. It has been well investigated in the study. For many details, see monographs, for example, [1, 3] and many others.

Let us consider that velocities in constraints Eq. (2b) are increased by virtual increments $\delta \dot{\mathbf{x}}$, so that they read

$$\sum _{s=1}^n{A}_{\mathit{rs}}\left({\dot{x}}_s+\delta {\dot{x}}_s\right)+B_r=0,r=1,2,..,l,$$E6

Deducting from Eq. (6), the initial state (Eq. (2b)) holds

$$\sum _{s=1}^n{A}_{\mathit{rs}}\delta {\dot{x}}_s=0,r=1,2,..,l.$$E7

Virtual increments of velocities $\delta \dot{\mathbf{x}}$ fit into constraints requested for constraints (Eq. (4)) and, consequently, in Eq. (5) the $\delta x_s$ can be replaced by $\delta {\dot{x}}_s$:

$$\sum _{s=1}^n\left(m_s{\ddot{x}}_s-X_s\right)\delta {\dot{x}}_s=0.$$E8

We revisit Eq. (2b) and perform differentiation with respect to $t$:

$$\sum _{s=1}^n\left(A_{\mathit{rs}}{\ddot{x}}_s+\frac{\mathrm{d}{A}_{\mathit{rs}}}{\mathrm{d}t}{\dot{x}}_s\right)+\frac{\mathrm{d}{B}_r}{\mathrm{d}t}=0,r=1,..,l,$$E9

where $\mathrm{d}/\mathrm{d}t$ represents the operator $\partial /\partial t+\sum _{i=1}^n{\dot{x}}_i\partial /\partial x_i$. Considering two possible movements of the system in identical initial state and velocities in time $t$, but with different accelerations $\ddot{\mathbf{x}}$ and $\ddot{\mathbf{x}}+\delta \ddot{\mathbf{x}}$, then we obtain with respect to Eq. (9):

$$\sum _{s=1}^n{A}_{\mathit{rs}}\delta {\ddot{x}}_r=0,r=1,..,l.$$E10

It means that virtual accelerations $\delta \ddot{\mathbf{x}}$ satisfy constraints requested similarly like virtual displacements or velocities following Eqs. (4) or (7). Therefore, we can write:

$$\sum _{s=1}^n\left(m_r{\ddot{x}}_s-X_s\right)\delta {\ddot{x}}_s=0.$$E11

Some authors call relations (Eqs. (5), (8), and (11)) as the first, second, and third form of the system equation, see, for example, [3] and others.

Let us multiply each equation in the system (Eq. (2a)) by velocities ${\dot{x}}_s$. Summing them together, one obtains

$$\sum _{s=1}^n{m}_s{\ddot{x}}_s{\dot{x}}_s=\sum _{s=1}^n{X}_s{\dot{x}}_s+\sum _{r=1}^l\sum _{s=1}^n{\lambda }_r{A}_{\mathit{rs}}{\dot{x}}_s,$$E12

which can be rewritten in the form

$$\frac{\mathrm{d}\mathcal{T}}{\mathrm{d}t}+\frac{\mathrm{d}\mathcal{V}}{\mathrm{d}t}=\sum _{s=1}^n{\tilde{X}}_s{\dot{x}}_s,$$E13

where $\mathcal{T},\mathcal{V}$ are kinetic and potential energy, respectively, and ${\tilde{X}}_s$ are forces, which cannot be included into the potential energy $\mathcal{V}$. In other words, relation Eq. (13) indicates that the change of the full energy (kinetic and potential) is equivalent to power (work on velocities) of all forces ${\tilde{X}}_s$, which do not contribute to the potential energy $\mathcal{V}$. Relation Eq. (13) corresponds to equilibrium condition Eq. (1), where functions of excitation and dissipation $\mathcal{P}\left(t\right),\mathcal{S}\left(t\right)$ correspond with the influence of non-potential forces ${\tilde{X}}_s$.

3. Lagrange’s equations

The original coordinates $\mathbf{x}$ should be replaced with respect to Lagrangian coordinates $\mathbf{q}$. The reason is that they represent the most inherent coordinates respecting the real movement of the system and configuration of external forces. Let us write basic coordinates as functions of Lagrangian ones:

$$x_r=x_r\left(q_1,..,q_n,t\right),r=1,..,n.$$E14

It can be easily shown

$$\frac{\partial {\dot{x}}_r}{\partial {\dot{q}}_s}=\frac{\partial x_r}{\partial q_s},\frac{\partial {\dot{x}}_r}{\partial q_s}=\frac{\mathrm{d}}{\mathrm{d}t}\left(\frac{\partial x_r}{\partial q_s}\right).$$E15

We reconsider Eq. (5) where the virtual displacement $\delta \mathbf{x}$ is replaced by

$$\delta x_r=\sum _{s=1}^n\frac{\partial x_r}{\partial q_s}\delta q_s,r=1,..,n\Rightarrow \sum _{s=1}^n\left\{\sum _{r=1}^n\left(m_r{\ddot{x}}_r-X_r\right)\frac{\partial x_r}{\partial q_s}\right\}\delta q_s=0.$$E16

The last Eq. (16) can be modified using Eqs. (15), which implies

$$\sum _{s=1}^n\left\{\sum _{r=1}^n{m}_r\left[\frac{\mathrm{d}}{\mathrm{d}t}\left({\dot{x}}_r\frac{\partial {\dot{x}}_r}{\partial {\dot{q}}_s}\right)-{\dot{x}}_r\frac{\partial {\dot{x}}_r}{\partial q_s}\right]-\sum _{r=1}^n{X}_r\frac{\partial x_r}{\partial q_s}\right\}\delta q_s=0,$$E17

This equation can be rewritten now in the form:

$$\sum _{s=1}^n\left\{\frac{\mathrm{d}}{\mathrm{d}t}\left(\frac{\partial \mathcal{L}}{\partial {\dot{q}}_s}\right)-\frac{\partial \mathcal{L}}{\partial q_s}-Q_s\right\}\delta q_s=0,$$E18

where it has been denoted:

$$\mathcal{L}=\sum _{r=1}^n{\dot{x}}_r\frac{\partial {\dot{x}}_r}{\partial q_s},Q_s=\sum _{r=1}^n{X}_r\frac{\partial x_r}{\partial q_s}.$$E19

Inspecting the polynomial $\mathcal{L}$, we recognize that it consists of the polynomial of the second and first degrees of components $\dot{\mathbf{q}}$ (coefficients are still functions of displacements $\mathbf{q}$ and time $t$) and the absolute part without any velocity components $\dot{\mathbf{q}}$. We can now assign the first part to the kinetic energy $\mathcal{T}$, while the part without velocities to the potential energy $\mathcal{V}$. So that $\mathcal{L}$ can be understood as the Lagrange function as usually defined

$$\mathcal{L}=\mathcal{T}-\mathcal{V},$$E20

provided the dynamic system studied is holonomic and no constraints are applied. In such a case, all variations $\delta q_s$ are independent and Eq. (18) can be fulfilled only if every coefficient in curly brackets vanishes individually. Consequently, we obtain Lagrange’s equations in the form:

$$\frac{\mathrm{d}}{\mathrm{d}t}\left(\frac{\partial \mathcal{T}}{\partial {\dot{q}}_s}\right)-\frac{\partial \mathcal{T}}{\partial q_s}+\frac{\partial \mathcal{V}}{\partial q_s}=Q_s,s_1,..,n,$$E21

where $Q_s$ are generalized external forces as functions of $\mathbf{q}$ and $t$. These forces are basically linear transforms of original forces $X_r$, see Eq. (19). If holonomic constraints are inserted, then the number of remaining degrees of freedom is lower $\left(k=n-l\right)$. Nevertheless, if there is possibility to define the system after elimination of inactive DOFs, then we can consider formally $k=n$ again and Eq. (21) remains in force.

Let us suppose now that the system includes $l$ non-holonomic constraints and those holonomic, which cannot be eliminated. Whatever is the reason for that, it still holds $k+l=n$. These constraints are described by constraints in Lagrange’s coordinates (analogous with Eq. (2b)) as follows:

$$\sum _{s=1}^n{C}_{\mathit{rs}}{\dot{q}}_s+D_r=0,r=1,..,l,\mathbf{D}=\left[D_r\right],\mathbf{D}\in {\mathbb{R}}^l,\mathbf{C}=\left[C_{\mathit{rs}}\right],\mathbf{C}\in {\mathbb{R}}^{l\times n},$$E22

This time, the variations $\delta q_s$ are not fully independent and only those components, which satisfy conditions:

$$\sum _{s=1}^n{C}_{\mathit{rs}}\delta q_s=0r=1,2,..,l,$$E23

can be regarded as independent.

In such a case, the right side of Eq. (21) should be completed:

$$\frac{\mathrm{d}}{\mathrm{d}t}\left(\frac{\partial \mathcal{T}}{\partial {\dot{q}}_s}\right)-\frac{\partial \mathcal{T}}{\partial q_s}+\frac{\partial \mathcal{V}}{\partial q_s}=Q_s+\sum _{r=1}^l{\lambda }_r{C}_{\mathit{rs}},s=1,..,n.$$E24

To the system, Eq. (24) should be attached $l$ constraints Eq. (22). So that, finally we have the system of $n+l$ equations with unknowns $\mathbf{q}$ and $\mathit{\lambda }$. Multipliers $\mathit{\lambda }$ are linearly related with forces in constraints. In particular cases, multipliers $\mathit{\lambda }$ can be physically interpreted, for instance, they can have a meaning of reactions of a body moving along a given trajectory. Very knowledgeable explanation about manipulation and interpretation of Lagrange’s multipliers from the viewpoint of a general theory as well as of employment in particular cases can be found in the monograph concerning non-holonomic systems, see [2]. For additional information and a large overview of additional literature resources, see [5].

The real dynamic system is always influenced by energy dissipation. Some simple models can be introduced using Rayleigh function $\mathcal{R}$, see, for example, [3]. This way is typically applicable, if linear viscous damping is considered and the Rayleigh function has a quadratic form in velocities ${\dot{q}}_s$. We can include this factor symbolically into Eq. (24), which reads now:

$$\frac{\mathrm{d}}{\mathrm{d}t}\left(\frac{\partial \mathcal{T}}{\partial {\dot{q}}_s}\right)-\frac{\partial \mathcal{T}}{\partial q_s}+\frac{\partial \mathcal{V}}{\partial q_s}=Q_s+\frac{\mathrm{\partial ℛ}}{\partial {\dot{q}}_s}+\sum _{r=1}^l{\lambda }_r{C}_{\mathit{rs}},s_1,..,n.$$E25

Hence, the completed system Eqs. (22) and (25) with $n+l$ unknowns can be considered. However, we should be aware that this supplement is rather intuitive and does not follow from any rigorous derivation, although in practice it is widely and successfully used. Nevertheless, comparison of this system with the general relation Eq. (1) introduced in Section 1 is obvious. Take a note that more sophisticated versions of Lagrange’s equations have been developed inspired by physical problems; see, for instance, generalized Lagrange-d’Alembert-Poincaré equations discussed in [11].

Let us add that many details (internal mechanisms and inclusion into governing system) concerning more sophisticated models of damping can be found in monographs of the rational dynamics, for example, [1, 4]. See also papers oriented to practical aspects of the damping either of natural, for example, rheological, aeroelastic origin, or intentionally included in order to achieve the highest damping effectiveness, for instance [13].

4. Appell-Gibbs function and equation system

Although the Appell-Gibbs approach is not referred so often in the study as the Lagrangian procedure, there are some monographs treating the analytical dynamics, for example, [1, 3], where detailed features of this method are explained. Moreover, journal papers can be found where special aspects of the Appell-Gibbs approach are discussed. A close relation of the fifth Gaussian form and the Gibbs equations from the viewpoint of Dynamics is studied, for example, [14, 15], important remarks for application are concerned in [16, 17], as well as possibilities of extension for systems with time-dependent masses [18] are indicated.

Let us briefly outline principal steps leading to the Appell-Gibbs differential system with respect to essentials ascertained and introduced in Section 2. We should be aware that generalized external forces $Q_s$, introduced in Eq. (19), follow in principle only $k$ degrees of freedom, which remained free; thereafter, $l$ constraints have been applied and the original number $n$ of DOFs has been reduced to $k=n-l,0<l\le n$. However, due to complicated relations inside the dynamic system, this fact is rather impossible to be employed in basic coordinates $x_s,s=1,..,n$ and Lagrange’s coordinates $q_s,s=1,..,n$ should be addressed, as we have also seen in previous Section 3. Nevertheless, it is worthy to involve only such coordinates $q_s$, which correspond to $k$ remained DOFs. It can be easily expressed in Lagrange’s coordinates, unlike basic coordinates $x_s$. So that, as the first step, we reformulate some expressions of Section 2 concerning the transform from basic to Lagrange’s coordinates.

Velocities ${\dot{x}}_r,r=1,..,n$ should be evaluated with respect to the fact that coordinates $x_r$ are functions of all Lagrange’s coordinates $q_s,s=1,..,k$ and time $t$, see Eq. (14):

$${\dot{x}}_r=\sum _{s=1}^k{\alpha }_{\mathit{rs}}{\dot{q}}_s+{\alpha }_r,r=1,..,n,\text{where}:{\alpha }_{\mathit{rs}}=\frac{\partial x_r}{\partial q_s},{\alpha }_r=\frac{\partial x_r}{\partial t},$$E26

which also implies

$$\delta {\dot{x}}_r=\sum _{s=1}^k{\alpha }_{\mathit{rs}}\delta {\dot{q}}_s,r=1,..,n,$$E27

Differentiation of Eq. (26) with respect to time gives

$${\ddot{x}}_r=\sum _{s=1}^k{\alpha }_{\mathit{rs}}{\ddot{q}}_s+\sum _{s=1}^k\frac{\mathrm{d}{\alpha }_{\mathit{rs}}}{\mathrm{d}t}{\dot{q}}_s+\frac{\mathrm{d}{\alpha }_r}{\mathrm{d}t},\frac{\mathrm{d}}{\mathrm{d}t}=\frac{\partial }{\partial t}+\sum _{m=1}^k{\dot{q}}_m\frac{\partial }{\mathrm{d}t},r=1,..,n.$$E28

The incremented acceleration vector, when keeping velocities and displacements, can be formulated as follows:

$${\ddot{x}}_r+\delta {\ddot{x}}_r=\sum _{s=1}^k{\alpha }_{\mathit{rs}}\left({\ddot{q}}_s+\delta {\ddot{q}}_s\right)+\sum _{s=1}^k\frac{\mathrm{d}{\alpha }_{\mathit{rs}}}{\mathrm{d}t}{\dot{q}}_s+\frac{\mathrm{d}{\alpha }_r}{\mathrm{d}t},r=1,..,n.$$E29

Deducting Eq. (28) from Eq. (29), one obtains

$$\delta {\ddot{x}}_r=\sum _{s=1}^k{\alpha }_{\mathit{rs}}\delta {\ddot{q}}_s,r=1,..,n.$$E30

Hence, it can be written

$$\sum _{r=1}^n{X}_r\delta {\ddot{x}}_r=\sum _{s=1}^k\left(\sum _{r=1}^n{X}_r{\alpha }_{\mathit{rs}}\right)\delta {\ddot{q}}_s=\sum _{s=1}^k{Q}_s\delta {\ddot{q}}_s,$$E31

where $Q_s$ are identical generalized forces, as they have been defined in Eq. (19). With reference to Eq. (11), we can reformulate this equation as follows:

$$\sum _{r=1}^n{m}_r{\ddot{x}}_r\delta {\ddot{x}}_r-\sum _{s=1}^k{Q}_s\delta {\ddot{q}}_s=0.$$E32

This relation will be used later, see Eq. (36).

As a principal step of this section, we define now the Gibbs function $\mathcal{G}$ concentrating “acceleration energy” included in all $n$ DOFs as follows:

$$\mathcal{G}=\frac{1}{2}\sum _{r=1}^n{m}_r{\ddot{x}}_r^2,$$E33

When we pass from basic to Lagrange’s coordinates, only $k$ active coordinates remain in force and so the expression Eq. (33) can be rewritten:

$$\mathcal{G}=\frac{1}{2}\sum _{r=1}^k{m}_r{\ddot{q}}_r^2,$$E34

Expressions Eqs. (33) and (34) differ only in terms independent from accelerations.

Let us introduce the function $\mathcal{H}$:

$$\mathcal{H}=\mathcal{G}-\sum _{s=1}^k{Q}_s{\ddot{q}}_s,$$E35

and evaluate its virtual increment:

$$\begin{array}{ll}\delta \mathcal{H}=\delta \left(\mathcal{G}-\sum _{s=1}^k{Q}_s{{\displaystyle \ddot{q}}}_s\right)& =\frac{1}{2}\sum _{r=1}^n{m}_r{\left({{\displaystyle \ddot{x}}}_r+\delta {{\displaystyle \ddot{x}}}_r\right)}^2-\frac{1}{2}\sum _{r=1}^n{m}_r{{\displaystyle \ddot{x}}}_r^2-\sum _{s=1}^k{Q}_s\delta {{\displaystyle \ddot{q}}}_s\\ & =\frac{1}{2}\sum _{r=1}^n{m}_r{\left(\delta {{\displaystyle \ddot{x}}}_r\right)}^2+\left(\sum _{r=1}^n{m}_r{{\displaystyle \ddot{x}}}_r\delta {{\displaystyle \ddot{x}}}_r-\sum _{s=1}^k{Q}_s\delta {{\displaystyle \ddot{q}}}_s\right).\end{array}$$E36

The last parenthesis vanishes due to the relation Eq. (32). Therefore, if $\delta \ddot{\mathbf{x}}\ne 0$, then the function $\delta \mathcal{H}$ is always positive:

$$\delta \mathcal{H}=\delta \left(\mathcal{G}-\sum _{s=1}^k{Q}_s{\ddot{q}}_s\right)>0,$$E37

which implies that accelerations ${\ddot{q}}_s,s=1,..,k$ should lead to a minimum of the function $\mathcal{H}$, which means:

$$\frac{\partial \mathcal{G}}{\partial {\ddot{q}}_r}=Q_r,r=1,..,k.$$E38

The energy dissipation terms $R_x,R_y,R_z$ should be added to the right side of Eq. (38). At this moment, the conformity of Eq. (38) with the equivalence Eq. (1) is well pronounced, similar like in the previous section. The system Eq. (38) should be completed by geometric constraints:

$${\dot{q}}_r=\sum _{s=1}^k{\beta }_{\mathit{rs}}{\dot{q}}_s+{\beta }_{r}r=k+1,..,n.$$E39

Equations (38) and (39) are the Gibbs-Appell differential system including $n$ equations, which can be written in the normal form and hence it is suitable to be immediately investigated using common methods.

The differential system (Eqs. (38) and (39)) represents the simplest and in the same time the most general form of equations of the dynamic system movement. The form of this system is very simple, and it can be used with the same effectiveness to the investigation of holonomic as well as non-holonomic systems, as the constraints can represent non-holonomic but also holonomic type of constraints. Unlike the Lagrangian approach, the non-holonomic or non-eliminable constraints do not augment the number of differential equations.

Procedure of the Appell-Gibbs equations employment in particular cases is obvious, looking back at this section. In the first step, the “so called kinetic energy of accelerations” $\frac{1}{2}{\sum }_{r=1}^N{m}_r{\ddot{x}}_r^2$ is composed using $n$ acceleration components of the vector $\ddot{\mathbf{x}}$. It represents the Appell-Gibbs function $\mathcal{G}$. In a general case, this function includes also all coordinates $\mathbf{x}$ and velocities $\dot{\mathbf{x}}$. Nevertheless, it is important that $\mathcal{G}$ in Lagrange’s coordinates contains only $k$ selected components of accelerations $\ddot{\mathbf{q}}$. Anyway, all $n$ components of $\dot{\mathbf{q}}$ and $\mathbf{q}$ are still included as a result of a transformation from basic to Lagrange’s coordinates.

It is worthy to remind that the differentiation outlined in Eq. (38) is very easy in a particular case. Indeed, let us realize that $\mathcal{G}$ can be symbolically expressed as a sum of quadratic function of accelerations ${\ddot{q}}_s,s=1,..,k\to {\mathcal{G}}_2$, linear function of these components ${\mathcal{G}}_1$ and function without accelerations ${\mathcal{G}}_0$. Differentiating ${\mathcal{G}}_2$, one obtains the relevant acceleration component in a linear form, which will be moved onto the left side together with a coefficient, which can be a function of all velocities and displacements ${\dot{q}}_s,q_s,s=1,..,n$. Differentiation of ${\mathcal{G}}_1$ leads to acceleration-free coefficients and ${\mathcal{G}}_0$ can be omitted leading to zeroes. Sometimes, the so-called reduced Appell-Gibbs function ${\mathcal{G}}^{\ast }$ is defined where ${\mathcal{G}}_0$ is a priori omitted.

In the second step, the work of $k$ given forces $\mathbf{Q}$ on $k$ virtual displacements $\mathbf{q}$ is carried out. It has the form $\sum _{s=1}^k{Q}_s\delta q_s$. We substitute now back into Eqs. (38) and add $l=n-k$ geometric constraints following Eqs. (39). So we obtain $k+l=n$ differential equations for $n$ components of the vector $\ddot{\mathbf{q}}\left(t\right)$. Take a note that no unknown multipliers $\mathit{\lambda }$ emerge here, which on the other hand increases the number of unknowns in a Lagrangian approach.

The procedure working with accelerations instead with velocities provides much simpler governing differential system. Unlike velocities, the acceleration components in the Appell-Gibbs function are included only in a few parts of energy expression. Therefore, all parts including only velocity and displacement components disappear during the differentiation of the Appell-Gibbs function with respect to ${\ddot{q}}_r,r=1,..,k$, and therefore they can be considered beforehand as unimportant.

Investigating problems with rotations, we work with Lagrange’s coordinates $\mathit{\omega }$, which represent in fact velocities. So that by solving the abovementioned differential system, the displacements and velocities $\mathit{\omega }$ emerge as results. Rotations themselves remain unattended. May be, it is a forfeit for a relative simplicity of the governing system in comparison with the Lagrangian approach. However, this shortcoming is mostly apparent only. The main part of the result represents usually displacement components, which are obtained without restrictions. Together with velocities $\mathit{\omega }$, they represent a full set of information needed to get through the shape of trajectories of the system response including rotation (illustrative example will be presented later in Section 6). If detailed rotations (not only velocities) are still needed, a subsequent integration can be performed independently using differential relations between rotation velocity vector $\mathit{\omega }$ and (for instance) Euler angles, see monographs [1, 3, 4] and others. They provide a detailed description of time history of a body orientation as a function of time $t$. This step can be useful, for instance, when a detailed animation is needed for presentation purposes.

5. Planar movement of a ball in a spherical cavity

5.1. Engineering motivation

Passive vibration absorbers of various types are very widely used in civil engineering. TV towers, masts, and other slender structures exposed to wind excitation are usually equipped by such devices. Conventional passive absorbers are of the pendulum type. Although they are very effective and reliable, they have several disadvantages limiting their application.

These shortcomings can be avoided using the absorber of the ball type. The basic principle comes out of a rolling movement of a metallic ball of a radius $r$ inside of a rubber-coated cavity of a radius $R>r$. This system is closed in an airtight case, see, for instance, Figure 1. First papers dealing with the theory and practical aspects of ball absorbers have been published during the last decade, see [13, 19].

Figure 1.

Dynamic scheme of (a) spherical pendulum absorber, (b) ball absorber, and (c) ball absorber during testing in a dynamic laboratory, see [19].

5.2. Planar layout of the system, Lagrangian procedure

The version, when the ball is forced to move solely in a vertical plane, has been thoroughly studied using Lagrangian approach in [20, 21] and other detailed papers dealing not only with theoretical aspects but also with experimental verification in the laboratory and in situ examining absorbers installed on real structures.

The cavity is fixed to a vibrating structure. Their dynamic character is represented by a linear single degree of freedom (SDOF) system represented by a mass $M$. Inside of the cavity, the ball $m$ in a vertical plane is moving, that is two degrees of freedom (TDOF) system should be investigated, as it is outlined in Figure 2. It follows from geometric relations:

$$R\cdot \phi =r\left(\psi +\phi \right)\Rightarrow \mathit{r\psi }={\varrho }_r\phi ,$$E40

where ${\rho }_r=R-r$. It holds for vertical or horizontal components of a displacement and velocity of the internal ball center:

$$\left.\begin{array}{rrr}\text{horiz}.:& u+{\mathrm{\varrho }}_r\cdot \text{sin}\phi \Rightarrow & \dot{u}+{\mathrm{\varrho }}_r\dot{\phi }\text{cos}\phi ,\\ \text{vert}.:& {\mathrm{\varrho }}_r\cdot \text{cos}\phi \Rightarrow & -{\mathrm{\varrho }}_r\dot{\phi }\text{sin}\phi .\end{array}\right\}$$E41

Figure 2.

Basic scheme of a system.

Kinetic energy of a moving system of the ball $m$ and the cavity $M$ can be written in a form:

$$\begin{array}{c}\mathcal{T}=\frac{1}{2}m\left[{\left(\dot{u}+{\mathrm{\varrho }}_r\dot{\phi }\text{cos}\phi \right)}^2+{\mathrm{\varrho }}_r^2{\dot{\phi }}^2{\text{sin}}^2\phi \right]+\frac{1}{2}J{\dot{\psi }}^2+\frac{1}{2}M{\dot{u}}^2=\frac{1}{2}\left(m+M\right){\dot{u}}^2+m{\mathrm{\varrho }}_r\dot{u}\dot{\phi }\text{cos}\phi \\ +\frac{m}{2\kappa }{\mathrm{\varrho }}_r^2{\dot{\phi }}^2,\end{array}$$E42

where $m/\kappa =m+J/r^2\Rightarrow \kappa =5/7$, while the potential energy is given by an expression:

$$\mathcal{V}=\mathit{mg}{\mathrm{\varrho }}_r\left(1-\text{cos}\phi \right)+\frac{1}{2}{\mathit{Cu}}^2.$$E43

The damping should be introduced in a form of a simple Rayleigh function:

$$\mathcal{R}=\frac{1}{2}\left({\mathit{Mb}}_u{\dot{u}}^2+{\mathit{mb}}_{\phi }{\mathrm{\varrho }}_r^2{\dot{\phi }}^2\right).$$E44

$m,M$ – mass of the ball $m$, mass of the cavity, $M$ representing the protected structure;

$J$ – inertia moment of the ball $m$;

$b_u,b_{\phi }$ – damping coefficients (logarithmic decrements, linear viscous damping);

Expressions Eqs. (42), (43), and (44) should be put into Lagrange’s equations of the second type, see Eqs. (24) or (25) and monographs, for example, [1, 3, 4] and others:

$$\begin{array}{l}\sum _{r=1}^n\left\{\frac{d}{\mathit{dt}}\left(\frac{\partial \mathcal{T}}{\partial {\dot{q}}_r}\right)-\frac{\partial \mathcal{T}}{\partial q_r}+\frac{\partial \mathcal{V}}{\partial q_r}+\frac{\mathrm{\partial ℛ}}{\partial {\dot{q}}_r}\right\}\delta q_r=P_r\left(t\right),\\ q_1=u=\zeta \cdot {\mathrm{\varrho }}_r,q_2=\phi ,P_u\left(t\right)=p\left(t\right)\cdot M{\mathrm{\varrho }}_r,P_{\phi }\left(t\right)=0,\end{array}$$E45

which give the governing equations of the system:

$$\begin{array}{rrrr}\ddot{\phi }+\kappa b_{\phi }\dot{\phi }+{\mathit{\kappa \omega }}_m^2\text{sin}\phi +\kappa \ddot{\zeta }\cdot \text{cos}\phi & =& 0,& \left(a\right)\\ \mu \ddot{\phi }\text{cos}\phi -\mu {\dot{\phi }}^2\text{sin}\phi +\left(1+\mu \right)\ddot{\zeta }+b_u\dot{\zeta }+{\omega }_M^2\zeta & =& p\left(t\right),& \left(b\right)\\ \mu =m/M,{\omega }_M^2=C/M,{\omega }_m^2=g/{\mathrm{\varrho }}_r.& & & \left(c\right)\end{array}$$E46

Equation (46) describes 2D movement of a ball absorber under excitation by the force $P\left(t\right)$ at any arbitrary deviation amplitudes including incidental transition through a limit cycle toward an open regime.

5.3. Illustration of some planar system features

Analysis of the governing system (Eqs. (46)) has been done in a couple of papers, for example, [20, 21]. Investigation has been carried out using the harmonic or multi-harmonic balance method, see, for example, [22, 23], respectively.

The system is auto-parametric, see, for example, [24] and other resources. Very rich overview of a theoretical basis of auto-parametric systems can be found in [25]. Expecting a single mode response, the Harmonic balance-based methods are applicable. Following approximate expressions for excitation and response can be written (cf., e.g., [22]):

$$\begin{array}{ll}p\left(t\right)& =p_0\text{sin}\left(\mathit{\omega t}\right),\\ \phi \left(t\right)& =\alpha \text{sin}\left(\mathit{\omega t}\right)+\beta \text{cos}\left(\mathit{\omega t}\right),\\ \zeta \left(t\right)& =\gamma \text{sin}\left(\mathit{\omega t}\right)+\delta \text{cos}\left(\mathit{\omega t}\right).\}\end{array}$$E47

Having four new variables $\alpha =\alpha \left(t\right),\beta =\beta \left(t\right),\gamma =\gamma \left(t\right),\delta =\delta \left(t\right)$ instead of two original unknowns $\phi \left(t\right),\zeta \left(t\right)$, two additional conditions can be freely chosen:

$$\dot{\alpha }\text{sin}\left(\mathit{\omega t}\right)+\dot{\beta }\text{cos}\left(\mathit{\omega t}\right)=0,\dot{\gamma }\text{sin}\left(\mathit{\omega t}\right)+\dot{\delta }\text{cos}\left(\mathit{\omega t}\right)=0.$$E48

After substituting Eqs. (47) and (48) into Eqs. (46) and substituting the $\text{sin}\varpi$ and $\text{cos}\phi$ functions by two terms of Taylor expansion, the harmonic balance procedure gives the differential system for unknown amplitudes $\mathbf{Z}={\left(\alpha ,\beta ,\gamma ,\delta \right)}^T$, see, for example, [21, 23]:

$$\mathbf{M}\left(\mathbf{Z}\right)\dot{\mathbf{Z}}=\mathbf{F}\left(\mathbf{Z}\right).$$E49

System (49) for amplitudes $\mathbf{Z}\left(t\right)$ is meaningful if they are functions of a “slow time,” in other words, if their changes within one period $2\pi /\omega$ are small or vanishing and individual steps of the harmonic balance operation are acceptable. The matrix $\mathbf{M}$ and the right-hand side vector $\mathbf{F}$ have the following form:

$$\mathbf{M}=\left(\begin{array}{llll}0& -\omega & -\frac{1}{4}\mathit{\text{αβκω}}& \frac{1}{8}\mathit{\kappa \omega }{A}_{\alpha }\\ \omega & 0& -\frac{1}{8}\mathit{\kappa \omega }{A}_{\beta }& \frac{1}{4}\mathit{\text{αβκω}}\\ -\frac{1}{8}\mathit{\mu \omega }{A}_{\beta }& \frac{1}{4}\mathit{\text{αβμω}}& \left(\mu +1\right)\omega & 0\\ -\frac{1}{4}\mathit{\text{αβμω}}& \frac{1}{8}\mathit{\mu \omega }{A}_{\alpha }& 0& -\left(\mu +1\right)\omega \end{array}\right),$$E50

$$\mathbf{F}=\frac{1}{48}\left(\begin{array}{l}6A_0\kappa \left(3{\mathit{\gamma \omega }}^2-{\mathit{\alpha \omega }}_m^2\right)+12{\omega }^2\left(\kappa \left(\mathit{\alpha \beta \delta }+\left(8-{\beta }^2\right)\gamma \right)-4\alpha \right)-48\mathit{\beta \kappa \omega }{b}_{\phi }\\ 6A_0\kappa \left({\mathit{\delta \omega }}^2-{\mathit{\beta \omega }}_m^2\right)+12{\omega }^2\left(\mathit{\alpha \gamma \kappa }+\mathit{\beta \delta \kappa }-4\right)\beta +48\mathit{\alpha \kappa \omega }{b}_{\phi }\\ {\omega }^2\left(A_0\left(A_0+22\right)\mathit{\beta \mu }-16\left(3\delta \left(\mu +1\right)-4\mathit{\beta \mu }\right)\right)+48\left(\mathit{\gamma \omega }{b}_u+{\mathit{\delta \omega }}_M^2\right)\\ {\omega }^2\left(A_0\left(A_0+22\right)\mathit{\alpha \mu }-16\left(3\gamma \left(\mu +1\right)-4\mathit{\alpha \mu }\right)\right)-48\left(\mathit{\delta \omega }{b}_u-{\mathit{\gamma \omega }}_M^2+p_0\right)\end{array}\right),$$E51

where $A_0={\alpha }^2+{\beta }^2-8,A_{\alpha }=3{\alpha }^2+{\beta }^2-8,A_{\beta }={\alpha }^2+3{\beta }^2-8.$

Let us consider stationary response of the system. In this case, the derivatives $\mathrm{d}\mathbf{Z}/\mathrm{d}t$ vanish and the right-hand side has to vanish too. Eq. (49) degenerates to the form of

$$\mathbf{F}\left(\mathbf{Z}\right)=0$$E52

Thus, to identify the stationary solutions, the zero solution points of $\mathbf{F}$, depending on the excitation frequency and amplitude, should be traced. In the same time, the signum and the zero points of the Jacobian $\text{det}\left(\mathbf{JF}\right)$ have to be checked. The negative value of the Jacobian for a particular point indicates that the corresponding solution is stable, whereas when Jacobian vanishes a bifurcation could occur.

The curve $\mathbf{F}\left(\alpha ,\beta ,\gamma ,\delta ,\omega \right)=0$, projected into the planes $\left(\omega ,R\right)$ or $\left(\omega ,S\right)$ (for $S^2={\gamma }^2+{\delta }^2$), forms the resonance curves known from the analysis of linear oscillators. However, correspondence of this curve to the original Eq. (46) is limited to the case of stationary response. It is necessary to remind that limits of stationarity of the response cannot be determined from properties of Eq. (52) itself. The complete Eq. (49) has to be taken into account for this purpose.

With respect to actual experiences regarding passive vibration absorbers and some interesting properties of system (46), the following reference input data have been introduced:

$$M=10.0;m=2.0;{\mathrm{\varrho }}_r=0.71;b_{\phi }=0.1;b_u=0.2;C=140;p_o=0.5÷2.5.$$E53

Utilizing Eqs. (52) and (51), the nonlinear resonance curves describing the stationary response of system (46) can be obtained. A set of such curves for excitation amplitudes $p_0=\mathrm{0.25,0.5,1,1.5,2.5}$ is shown in Figure 3. It is obvious for the first view the nonlinear character manifesting oneself by a dependence of a position of extreme points on an amplitude of excitation force. This effect is visible predominantly in a neighborhood of a conventional “linear” natural frequency of the absorber, although also the second natural frequency corresponding to the original natural frequency of the structure is affected. The resonance curves are typical for a system with “softening” nonlinearities.

Figure 3.

Nonlinear resonance curves describing the stationary response of the system for excitation amplitudes $p_0=0.25,0.5,1,1.5,2.5$. Stable branches are shown as solid blue curves, unstable parts are indicated as the red dashed curves. Amplitudes, see Eq. (47), $R=\sqrt{{\alpha }^2+{\beta }^2}$ are shown in the left part of the figure, amplitudes $S=\sqrt{{\gamma }^2+{\delta }^2}$ are on the right.

6. Spatial version of the system, Appell-Gibbs procedure

6.1. Gibbs function

The spatial version of the ball absorber on the basis of rational dynamics has been widely investigated by authors of this chapter, see, for example, [26, 27]. Lagrangian approach and Appell-Gibbs procedures have been discussed in these papers combining analytical and numerical methods. Some important issues will be roughly outlined and for details see cited papers.

Unlike the planar version discussed in the previous section, the Appell-Gibbs approach is used to formulate the governing nonlinear differential system. The authors tried to formulate the spatial version using the Lagrangian procedure as well, see [28]. Although the governing system of the respective holonomic system has been successfully assembled, the further analysis appeared very cumbersome, and therefore, it has been given up to follow this way. Thus, the Appell-Gibbs approach is used to formulate the governing system. Its structure is much more transparent and represents a wider option of analytical-numerical investigation of detailed properties of the ball trajectories within the cavity.

With respect to Sections 2 and 4, the first step represents to construct the Appell-Gibbs function (often referred to as an energy acceleration function) defined as follows:

$$\mathcal{G}=\frac{1}{2}M\left({\ddot{u}}_{\mathit{Gx}}^2+{\ddot{u}}_{\mathit{Gy}}^2+{\ddot{u}}_{\mathit{Gz}}^2\right)+\frac{1}{2}J\left({\dot{\omega }}_x^2+{\dot{\omega }}_y^2+{\dot{\omega }}_z^2\right),$$E54

where $M$ is the mass of the ball, $J$ is central inertia moment of the ball with respect to point $G$, $\mathit{\omega }$ the angular velocity vector of the ball with respect to its center $G$, ${\mathbf{u}}_G$ the displacement of the ball center with respect to absolute origin $O$, $C$ contact point of the ball and cavity, $A$ moving origin related with the cavity in its bottom point, see Figure 4. Coordinates $\mathbf{x}=\left[x,y,z\right]$ are Cartesian coordinates with origin in the point $O$. Hence, it holds:

$$\begin{array}{ccc}{\mathbf{u}}_G={\mathbf{u}}_A+{\mathbf{u}}_C+{\mathbf{u}}_n,& {\mathbf{u}}_n=r\cdot \mathbf{n}& \\ {{\displaystyle \dot{\mathbf{u}}}}_G={{\displaystyle \dot{\mathbf{u}}}}_A+{\displaystyle \underbrace{{{\displaystyle \dot{\mathbf{u}}}}_C+r\cdot {{\displaystyle \dot{\mathbf{u}}}}_n,}}\\ \rho {{\displaystyle \dot{\mathbf{u}}}}_C,& \rho =1-r/R,& \left(\mathrm{cf}.\mathrm{Eq}.\left(40\right):{\mathrm{\varrho }}_r=R-r\right),\end{array}$$E55

where ${\mathbf{u}}_A$ is the displacement of the moving origin $A$ with respect to absolute origin $O$, ${\mathbf{u}}_C$ the displacement of the contact point $C$ with respect to moving origin $A$, ${\mathbf{u}}_n$ the displacement of the ball center $G$ with respect to contact point $C$, $\mathbf{n}$ the cavity normal unit vector in point $C$.

Figure 4.

Ball rotation vector in moving coordinates.

Geometry of the cavity (radius $R$) with respect to moving origin $A$ is given by equation:

$$x_A^2+y_A^2+{\left(z_A-R\right)}^2=R^2,$$E56

where ${\mathbf{x}}_A=\left[x_A,y_A,z_A\right]$ are Cartesian coordinates with origin in the moving origin $A$.

Using Pfaff theorem and adopting a conjecture of non-sliding contact between the ball and the cavity, the respective non-holonomic constraints of “perfect” rolling can be deduced after a longer manipulation:

$$\begin{array}{r}{\dot{u}}_{\mathit{Gx}}={\dot{u}}_{\mathit{Ax}}+\rho \left({\omega }_y\left(u_{\mathit{Cz}}-R\right)-{\omega }_z{u}_{\mathit{Cy}}\right),\\ {\dot{u}}_{\mathit{Gy}}={\dot{u}}_{\mathit{Ay}}+\rho \left({\omega }_z{u}_{\mathit{Cx}}-{\omega }_x\left(u_{\mathit{Cz}}-R\right)\right),\\ {\dot{u}}_{\mathit{Gz}}=+\rho \left({\omega }_x{u}_{\mathit{Cy}}-{\omega }_y{u}_{\mathit{Cx}}\right),\end{array}$$E57

where $\rho =1-r/R$.

In order to substitute for accelerations ${\mathbf{u}}_o$ into the Appell function (Eq. (54)), let us differentiate constraints Eqs. (57).

Several manipulations provide expressions for components of the ball center acceleration ${\ddot{\mathbf{u}}}_G$, which consist of acceleration in the moving origin $A:{\ddot{\mathbf{u}}}_A$ representing the given external kinematic excitation and acceleration related to the point $A$ being given by an expression: $\rho {\ddot{\mathbf{u}}}_C$:

$$\begin{array}{l}{\ddot{u}}_{\mathit{Gx}}={\ddot{u}}_{\mathit{Ax}}+\rho \left({\dot{\omega }}_y\left(u_{\mathit{Cz}}-R\right)-{\dot{\omega }}_z{u}_{\mathit{Cy}}\right)+\rho \left({\omega }_y{\dot{u}}_{\mathit{Cz}}-{\omega }_z{\dot{u}}_{\mathit{Cy}}\right),\\ {\ddot{u}}_{\mathit{Gy}}={\ddot{u}}_{\mathit{Ay}}+\rho \left({\dot{\omega }}_z{u}_{\mathit{Cx}}-{\dot{\omega }}_x\left(u_{\mathit{Cz}}-R\right)\right)+\rho \left({\omega }_y{\dot{u}}_{\mathit{Cz}}-{\omega }_z{\dot{u}}_{\mathit{Cy}}\right),\\ {\ddot{u}}_{\mathit{Gz}}={\ddot{u}}_{\mathit{Az}}+\rho \left({\dot{\omega }}_x{u}_{\mathit{Cy}}-{\dot{\omega }}_y{u}_{\mathit{Cx}}\right)+\rho \left({\omega }_x{\dot{u}}_{\mathit{Cy}}-{\omega }_y{\dot{u}}_{\mathit{Cx}}\right),\end{array}$$E58

Because the kinematic excitation is supposed to be horizontal, ${\ddot{u}}_{\mathit{Az}}=0$ into Eqs. (58) should be substituted.

Expressions Eqs. (58) are to be substituted into Eq. (54). Thereby, we obtain the Appell-Gibbs function $\mathcal{G}$ for the system investigated. The function $\mathcal{G}$ can be significantly simplified keeping only terms including second-time derivatives ${\ddot{\mathbf{u}}}_G$ and $\dot{\mathit{\omega }}$, which represent second-time derivatives of respective rotations. This step provides the reduced Appell-Gibbs function ${\mathcal{G}}^r$. Using ${\mathcal{G}}^r$, one can write the Appell-Gibbs differential system:

$$\partial {\mathcal{G}}^r/\partial {\dot{\omega }}_x=F_{\mathit{Gx}},\partial {\mathcal{G}}^r/\partial {\dot{\omega }}_y=F_{\mathit{Gy}},\partial {\mathcal{G}}^r/\partial {\dot{\omega }}_z=F_{\mathit{Gz}},$$E59

where ${\mathbf{F}}_G$ is the external force vector acting in ball center $G$. Vector ${\mathbf{F}}_G$ is determined subsequently using the virtual displacements principle. Let us introduce the quasi-coordinates ${\phi }_x,{\phi }_y,{\phi }_z$ where ${\omega }_x={\dot{\phi }}_x,{\omega }_y={\dot{\phi }}_y,{\omega }_z={\dot{\phi }}_z$. The only external force acting in the ball center is the gravity. Therefore, the elementary work performed can be expressed as

$$\delta F_G=-\mathit{mg}\cdot \delta u_{\mathit{Gz}}.$$E60

Virtual displacement $\delta u_{\mathit{Gz}}$ can be determined using the third non-holonomic constraint in Eqs. (57). It holds

$$\delta u_{\mathit{Gz}}=\rho \left(u_{\mathit{Cy}}{\mathit{\delta \phi }}_x-u_{\mathit{Cx}}{\mathit{\delta \phi }}_y\right),$$E61

and therefore

$$\delta F_G=-\mathit{mg\rho }\left(u_{\mathit{Cy}}{\mathit{\delta \phi }}_x-u_{\mathit{Cx}}{\mathit{\delta \phi }}_y\right).$$E62

At the same time, the elementary work can be expressed in terms of quasi-coordinates:

$$\delta F_G=F_{\mathit{Gx}}{\mathit{\delta \phi }}_x+F_{\mathit{Gy}}{\mathit{\delta \phi }}_y+F_{\mathit{Gz}}{\mathit{\delta \phi }}_z.$$E63

Comparing coefficients at respective virtual components ${\mathit{\delta \phi }}_x,{\mathit{\delta \phi }}_y,{\mathit{\delta \phi }}_z$, in Eqs. (61) and (63), one obtains

$$F_{\mathit{Gx}}=-\mathit{\rho mg}\cdot u_{\mathit{Cy}},F_{\mathit{Gy}}=\mathit{\rho mg}\cdot u_{\mathit{Cx}},F_{\mathit{Gz}}=0.$$E64

The damping will be introduced later in Section 6.3 in order to separate energy conservative approach and enable to discuss various stationary regimes with respect to parameter and excitation settings.

6.2. Governing system

Carrying out the differentiation outlined in Eqs. (59), and respecting Eqs. (64), it can be written after some adaptations:

$$J_s{{\displaystyle \dot{\omega }}}_x-u_{\mathit{Cx}}{{\displaystyle \dot{\mathrm{\Omega }}}}_s=\left(\left({{\displaystyle \ddot{u}}}_{\mathit{Ay}}+\rho \left({\omega }_z{{\displaystyle \dot{u}}}_{\mathit{Cx}}-{\omega }_x{{\displaystyle \dot{u}}}_{\mathit{Cz}}\right)\right)\left(u_{\mathit{Cz}}-R\right)\right.\left.-u_{\mathit{Cy}}\left(g+\rho \left({\omega }_x{{\displaystyle \dot{u}}}_{\mathit{Cy}}-{\omega }_y{{\displaystyle \dot{u}}}_{\mathit{Cx}}\right)\right)\right),$$

$$J_s{{\displaystyle \dot{\omega }}}_y-u_{\mathit{Cy}}{{\displaystyle \dot{\mathrm{\Omega }}}}_s=\left(-\left({{\displaystyle \ddot{u}}}_{\mathit{Ax}}+\rho \left({\omega }_y{{\displaystyle \dot{u}}}_{\mathit{Cz}}-{\omega }_z{{\displaystyle \dot{u}}}_{\mathit{Cy}}\right)\right)\left(u_{\mathit{Cz}}-R\right)\right.\left.+u_{\mathit{Cx}}\left(g+\rho \left({\omega }_x{{\displaystyle \dot{u}}}_{\mathit{Cy}}-{\omega }_y{{\displaystyle \dot{u}}}_{\mathit{Cx}}\right)\right)\right),$$E65

$$J_s{{\displaystyle \dot{\omega }}}_z-\left(u_{\mathit{Cz}}-R\right){{\displaystyle \dot{\mathrm{\Omega }}}}_s=\left(\left({{\displaystyle \ddot{u}}}_{\mathit{Ax}}+\rho \left({\omega }_y{{\displaystyle \dot{u}}}_{\mathit{Cz}}-{\omega }_z{{\displaystyle \dot{u}}}_{\mathit{Cy}}\right)\right)u_{\mathit{Cy}}\right.\left.-\left({{\displaystyle \ddot{u}}}_{\mathit{Ay}}+\rho \left({\omega }_z{{\displaystyle \dot{u}}}_{\mathit{Cx}}-{\omega }_x{{\displaystyle \dot{u}}}_{\mathit{Cz}}\right)\right)u_{\mathit{Cx}}\right),$$

where

$$\begin{array}{lll}{\dot{\mathrm{\Omega }}}_s& =& u_{\mathit{Cx}}{\dot{\omega }}_x+u_{\mathit{Cy}}{\dot{\omega }}_y+\left(u_{\mathit{Cz}}-R\right){\dot{\omega }}_z,\\ J_s& =& {\displaystyle \underbrace{\left(J+m{\rho }^2{R}^2\right)}}/m{\rho }^2.\\ & & \text{mass\hspace{0.17em}interia\hspace{0.17em}moment\hspace{0.17em}of\hspace{0.17em}the\hspace{0.17em}ball}\\ & & \text{with}\text{respect}\text{to}\text{center}\text{of}\text{the}\text{cavity}\end{array}$$E66

It can be shown that ${{\displaystyle \dot{\mathrm{\Omega }}}}_s=0$ and therefore the second column on the left side of the system Eqs. (65) should be omitted. External excitations are specified by movement or acceleration in the point $A$. Hence, kinematic excitation in $A$ is given as follows:

$$E_{\mathit{Ax}}\left(t\right)={\ddot{u}}_{\mathit{Ax}}/\rho ,E_{\mathit{Ay}}\left(t\right)={\ddot{u}}_{\mathit{Ay}}/\rho ,E_{\mathit{Az}}\left(t\right)=0,$$E67

as it can be seen in Eqs. (65). Provided we need to investigate the response processes in a vertical plane, only one component remains non-zero and the second vanishes as well.

In order to obtain the system Eqs. (65) in the form with first-time derivatives concentrated on the left side, the first derivatives ${\dot{\mathbf{u}}}_C$ in its right sides should be expressed in displacements ${\mathbf{u}}_C$ using non-holonomic constraints Eqs. (57):

$$\begin{array}{ll}{\dot{u}}_{\mathit{Cx}}& ={\omega }_y\left(u_{\mathit{Cz}}-R\right)-{\omega }_z{u}_{\mathit{Cy}},\\ {\dot{u}}_{\mathit{Cy}}& ={\omega }_z{u}_{\mathit{Cx}}-{\omega }_x\left(u_{\mathit{Cz}}-R\right),\\ {\dot{u}}_{\mathit{Cz}}& ={\omega }_x{u}_{\mathit{Cy}}-{\omega }_y{u}_{\mathit{Cx}}.\end{array}$$E68

Therefore, we obtained the system of six non-linear ODEs (Eqs. (65) and (68)) in a normal form with six unknown functions of time: $u_{\mathit{Cx}},u_{\mathit{Cy}},u_{\mathit{Cz}},{\omega }_x,{\omega }_y,{\omega }_z$. Vector ${\mathbf{u}}_C$ depicts displacements of the contact point and can be used to study the movement of the ball from a global point of view. Detailed behavior of the ball as a rotating body is given by angular velocities $\mathit{\omega }$. If the time history of rotation should be traced, then a subsequent run is necessary to obtain rotations by means of Euler angles as solution of the system of three ODEs with an input of angular velocities $\mathit{\omega }$.

6.3. Influence of the damping

Influence of the damping will be taken into account. Basically, two sources of the energy dissipation are ruling in the system: (1) dissipation due to air dynamic resistance and (2) energy loss in contact of the cavity and rolling ball. The former one can be neglected with respect to obvious geometric configuration of the device and relative velocity ball/cavity. Concerning the latter one, complicated energy dissipating processes are ruling in contact of the ball with cavity. Nevertheless, supposing that no slipping arises in the contact, the dissipation process can be approximated as proportional to relevant components of the angular velocity vector $\mathit{\omega }$ and the quality of the cavity/ball contact. Considering the obvious setting, the respective material coefficients characterizing the rolling movement of the ball can be considered constant regardless of the direction in the tangential plane to the cavity in the point $C$, see Figure 4. The coefficient determining the rotation resistance around the normal vector $\mathbf{n}$ in the contact point $C$ is different as a rule. Therefore, the resistance moment vector $\mathbf{D}$ can be expressed in moving coordinates $p,q,n$, see Figure 4, as follows:

$$\mathbf{D}={\left[D_p,D_q,D_n\right]}^T.$$E69

Components of the above vector can be written in a form as follows:

$$D_p={\kappa }_r\cdot {\omega }_p,D_q={\kappa }_r\cdot {\omega }_q,D_n={\kappa }_s\cdot {\omega }_n,$$E70

where ${\kappa }_r,{\kappa }_s$ are coefficients of “viscous resistance” of rolling and spinning. Their meaning is: the moment for a unity rotation per second, that is (Nms/rad).

Turning of the vector ${\mathbf{D}}_G={\left[D_{\mathit{Gx}},D_{\mathit{Gy}},D_{\mathit{Gz}}\right]}^T$ expressed in $\left(\mathit{xyz}\right)$ coordinates into the vector $\mathbf{D}$ can be written as

$$\mathbf{D}={\mathbf{T}}_C\cdot {\mathbf{D}}_G,$$E71

The transformation matrix $T_C$ reads

$${\mathbf{T}}_C=\left[\begin{array}{lll}\frac{x_C\left(-z_C+R\right)}{\mathit{R\nu }},& \frac{y_C\left(-z_C+R\right)}{\mathit{R\nu }},& \frac{\nu }{R}\\ \frac{-y_C}{\nu },& \frac{x_C}{\nu },& 0\\ \frac{-x_C}{R},& \frac{-y_C}{R},& \frac{-z_C+R}{R}\end{array}\right]$$E72

where ${\nu }^2=x_C^2+y_C^2$. The matrix ${\mathbf{T}}_C$ is orthogonal and, therefore, the inverse transformation goes using matrix ${\mathbf{T}}_C^{-1}={\mathbf{T}}_C^T$, in particular:

$${\mathbf{D}}_G={\mathbf{T}}_C^T\cdot \mathbf{D},$$E73

Components of the vector ${\mathbf{D}}_G$ should be incorporated onto the right side of Eqs. (59), where right sides should be completed. It means that the elementary work $\delta F_G$ following Eq. (60) must be completed by a negative dissipating work due to ${\mathbf{D}}_G$.

$$\delta F_G=-\mathit{mg}\cdot \delta u_{\mathit{Gz}}-{\mathbf{D}}_G\cdot \delta \mathit{\phi },$$E74

Repeating the further derivation like in Section 6.2, one can revisit the system Eqs. (65) and (68), where the right sides are completed and instead (Eqs. (65)) they read:

$$\begin{array}{r}{J}_s{\dot{\omega }}_x=(\left({\ddot{u}}_{\mathit{Ay}}+\rho \left({\omega }_z{\dot{u}}_{\mathit{Cx}}-{\omega }_x{\dot{u}}_{\mathit{Cz}}\right)\right)\left(u_{\mathit{Cz}}-R\right)-u_{\mathit{Cy}}\left(g+\rho \left({\omega }_x{\dot{u}}_{\mathit{Cy}}-{\omega }_y{\dot{u}}_{\mathit{Cx}}\right)\right))-D_{\mathit{Gx}}/m,\\ J_s{\dot{\omega }}_y=(-\left({\ddot{u}}_{\mathit{Ax}}+\rho \left({\omega }_y{\dot{u}}_{\mathit{Cz}}-{\omega }_z{\dot{u}}_{\mathit{Cy}}\right)\right)\left(u_{\mathit{Cz}}-R\right)+u_{\mathit{Cx}}\left(g+\rho \left({\omega }_x{\dot{u}}_{\mathit{Cy}}-{\omega }_y{\dot{u}}_{\mathit{Cx}}\right)\right))-D_{\mathit{Gy}}/m,\\ J_s{\dot{\omega }}_z=(\left({\ddot{u}}_{\mathit{Ax}}+\rho \left({\omega }_y{\dot{u}}_{\mathit{Cz}}-{\omega }_z{\dot{u}}_{\mathit{Cy}}\right)\right)u_{\mathit{Cy}}\left.-\left({\ddot{u}}_{\mathit{Ay}}+\rho \left({\omega }_z{\dot{u}}_{\mathit{Cx}}-{\omega }_x{\dot{u}}_{\mathit{Cz}}\right)\right)u_{\mathit{Cx}}\right)-D_{\mathit{Gz}}/m.\end{array}$$E75

Terms $D_{\mathit{Gx}}/m,D_{\mathit{Gy}}/m,D_{\mathit{Gz}}/m$ which are linear functions of ${\omega }_x,{\omega }_y,{\omega }_z$ determine the viscous type of the damping, although intensity in individual coordinates is variable depending on the position of the ball within the cavity.

6.4. Ball trajectories within the fixed cavity due to initial conditions

A large program of a ball trajectory investigation within a spherical cavity has been performed using the differential system (Eqs. (68) and (75)). Basically, it consists of two groups which are briefly illustrated in this and the next subsections. The first group concerns the fixed cavity (no excitation is applied). The only source of energy introduced is given by the initial deflection of the ball from equilibrium position in the point $A$ (“southern pole”), or in other words by non-homogeneous initial conditions.

Differential system (Eqs. (68) and (75)) admits a number of singular solutions which can serve as separating limits of zones within which regular solutions exhibit certain character of trajectory shape. Some of them can be found analytically from the differential system taking into account their special properties concerning individual response component along the trajectory as a whole or in certain points of these curves. For details, special papers should be referred. Take a note that most of them emerge when no damping is considered. The reason is that the trajectory should be quasi-periodic (or cyclic-stationary), which is impossible when damping is respected and no external energy supply is considered. Trajectories start in a certain point on a meridian into which the ball is elevated. Then, it is thrown horizontally along the cavity parallel circle. Let us mention a few of the most important:

1. circular trajectory in horizontal plane. No initial spin is considered $\left({\omega }_{n0}=0\right)$. The impulse applied corresponds with the initial velocity $\mathit{\omega }=\left[{\omega }_{\mathit{ps}}\mathrm{,0,0}\right]$, where it holds for ${\omega }_{\mathit{ps}}$:

$${\omega }_{\mathit{ps}}=\frac{\mathit{gM\rho }{u}_{\mathit{Cz}0}\left(2R-u_{\mathit{Cz}0}\right)}{\left(J+M{\rho }^2{R}^2\right)\left(R-u_{\mathit{Cz}0}\right)}.$$E76

This case is the most important and can be called separating circle (SC).

1. 2. circular trajectory in inclined plane, see Figure 5, ${\mathit{\omega }}_0=\left[\mathrm{100.0,0.0,0.0}\right]$. This state is exactly valid for ${\omega }_{p0}\to \infty$. The space spiral type trajectory changes from SC upwards successively into the upper hemisphere. Before the limit state for ${\omega }_{p0}\to \infty$ is reached, the osculation plane of the trajectory can be recognized. It rotates around the vertical axis with a descending angular velocity as far as it vanishes and osculation and operating planes coincide.

2. 3. trajectory of “kings crown form,” see Figure 5, ${\mathit{\omega }}_0=\left[\mathrm{5.817,0.0,5.0537}\right]$. Cases, when the initial spin is considered. For a special value of ${\omega }_{n0}=5.0537$ takes a shape visible in the picture. The apexes of this curve correspond to $\mathit{\omega }=\left[\mathrm{0.0,0.0},{\omega }_{n0}\right]$ and $\mathbf{u}=\left[\mathrm{0.0,0.0,0.0}\right]$, which is a clue to find forms and parameters of this special case. This trajectory is reached from SC, increasing the initial spin velocity until the limit value. If it is lower, the trajectory has the spiral form. For a higher value, it became a curly form, see Figure 5, ${\mathit{\omega }}_0=\left[\mathrm{5.817,0.0,10.0}\right]$. The limit state for infinite initial spin represents the ball apparently fixed in the initial point and not moving neither horizontally nor vertically.

Figure 5.

Illustration of the ball trajectories; cavity is not excited; energy only supply is due to non-homogeneous initial condition; in every triplet: movement time history of the contact point $C$: $u_{\mathit{Cx}},u_{\mathit{Cy}},u_{\mathit{Cz}}$; vertical view of trajectories $u_{\mathit{Cx}},u_{\mathit{Cy}}$ components; axonometric view of trajectories; parameters above triplets: initial values of ${\mathit{\omega }}_0=\left[{\omega }_{p0},{\omega }_{q0},{\omega }_{n0}\right]$ and damping parameters: ${\kappa }_r,{\kappa }_s$; line (a): no spin, no damping, line (b): spin considered, line (c): spin and damping considered.

Let us have a look at the bottom two pictures in Figure 5. They respect the influence of the damping. Coefficients ${\kappa }_r,{\kappa }_s$ are different as it corresponds to conditions in the real system. The left demonstrates trajectory for positive initial spin and the right for negative initial spin. The transition through limit cases mentioned earlier is visible. The trajectory obviously finishes in the bottom “southern pole” of the cavity.

6.5. Ball trajectories within kinematically excited cavity

The second group of tests deals with the cavity which is undergone to kinematic excitation in a horizontal plane (only one-direction excitation is reported here).

Two extensive series of tests demonstrate the auto-parametric character of the system. In the first series, the response has been evaluated separately for every excitation frequency $\omega$ starting from homogeneous initial conditions. Figure 6 shows some selected results of numerical simulations which follow from the differential system (Eqs. (68) and (75)). We briefly point out a couple of features visible in Figure 6. In the picture (a), we can see the maximal horizontal amplitude of the ball trajectory, when the cavity is kinematically excited in the horizontal plane in $x$ direction. The solid curve represents max $\mid u_{\mathit{Cx}}\mid$ and the dashed curve is max $\mid u_{\mathit{Cy}}\mid$ as functions of the exciting frequency $\omega$. We can see that in the interval $\omega \in \left(0,2.84\right)$, the semi-trivial solution is stable and so $u_y=0$. The point $\omega =2.84$ is a beginning of the resonance zone, which spans in $\omega \in \left(2.84,2.99\right)$, where auto-parametric resonance occurs and amplitudes of both response components are commeasurable. For $\omega >2.99$, the semi-trivial solution is regained. Samples of the trajectory shape are plotted in picture (b) for four frequencies $\omega =2.84,2.88,2.92,2.96$. Their vertical views demonstrate the character of the semi-trivial and the auto-parametric resonance states. Take a note that the trajectory since $\omega =2.94$ is a simple ellipse-like curve, which does not exhibit any symptom of a chaotic process. Compare this finding with analysis concerning the sweeping up and down excitation frequency for $\omega$ around and above $B_2$ bifurcation point (BP) (see Figure 7 and explanation later in this subsection).

Figure 6.

Response of the ball in the resonance and adjacent zones due to harmonic horizontal excitation of the cavity: (a) amplitude of the displacement as a function of the excitation frequency; (b) vertical views of the ball trajectory for frequencies $\omega =2.84,2.88,2.92,2.96$.

Figure 7.

Amplitudes of the ball displacement under cavity harmonic excitation, when the frequency is swept up and down: (a) amplitudes overview in the interval $\omega \in \left(1.0,8.0\right)$, (b) zooming in the interval $\omega \in \left(2.4,3.1\right)$; curves: solid red—max $\mid u_{\mathit{Cx}}\mid$, dashed red—max $\mid u_{\mathit{Cy}}\mid$, solid black—absolute displacement amplitude, blue dashed—attraction boundary between $b_{u1}$ and $b_{l1}$.

The second series has been controlled by sweeping the excitation frequency up and down in a large interval and in several detailed regimes in the area of the auto-parametric resonance zone. A few of the results are visible in Figure 7. Picture (a) demonstrates amplitudes max $\mid u_{\mathit{Cx}}\mid$ (solid curves) and max $\mid u_{\mathit{Cy}}\mid$ (dashed curve) and the total amplitude $u_{\mathit{Cr}}$ in the interval $\omega \in \left(1.0,8.5\right)$. Picture (b) is the magnified detail of picture (a) within the interval $\omega \in \left(2.80,3.05\right)$ in order to make visible the resonance zone.

Let us pay attention to bifurcation points (BPs). There are obviously concentrated in the resonance zone. In principle, they can be classified into two categories. The most important reveal $B_1$ and $B_2$. In the latter one, two branches start. The lower one $b_{l2}$ approaches zero for $\omega \to \infty$ which indicates the non-moving ball in the vertical view. This branch takes place in the vertical plane and basically has a form of semi-trivial solution. Its stability increases with rising $\omega >{\omega }_{B2}$ as it follows from decreasing negative values of the Lyapunov exponent and of inspection of the relevant stability basins. The upper branch $b_{u2}$ is spatial. It follows from the resonance zone where the spatial response type has a chaotic character. The relevant attractor reveals as an annular concentric area with diminishing width with increasing $\omega$. The trajectory very quickly approaches a circular form in the horizontal plane. Its level with respect to the vertical axis rises and approaches “equatorial” position. However, the stability of this trajectory decreases, and we can see in Figure 7 that around $\omega =8.0$ even numerical perturbations of the integration process can overcome the stability limit (despite very small integration step) and the response trajectory falls down to the lower branch in the point $D_2$. Its position is not fixed. If hypothetically zero perturbation occurs, it could shift to infinity and approach together with the branch $b_{u2}$ the asymptote at the level $R=1$. Observing black max $\mid u_{\mathit{Cx}}\mid$ and red max $\mid u_{\mathit{Cy}}\mid$ parts of $b_{u2}$, we can see that they are getting coincide with increasing $\omega$. It means that trajectory approaches the circle with radius $R=1$.

Let us briefly discuss the shape of the response amplitudes for $\omega$ below BP $B_1$ and $B_4$. The BP $B_1$ is reached sweeping up along the branch $b_{l1}$, when it loses planar character passing through $B_4$. In such a case, the spatial response type emerges, exhibiting a chaotic response since $B_1$. This fact is obvious also looking at the dashed red curve representing $u_{\mathit{Cy}}$, which is trivial as far as $B_4$ and can bring the system from the semi-trivial solution into the auto-parametric resonance starting $B_4$. Take a note that passing BP $B_4$, planar response can remain in force, if any perturbation is avoided. It meets in $B_3$ the branch $b_{l2}$ following also a planar path for swept up $\omega$. Branch $b_{u1}$ starts in $B_1$. Its stability rapidly decreases with descending $\omega$. The point $D_1$ illustrates its limited extent in sub-resonance zone. This feature is visible observing the curve $b_s$, which represents a limit separating the area of attraction to $b_{u1}$ and to $b_{l1}$. Take a note that $b_s$ starts in $B_4$ and approaches $D_1$, despite hypothetically it goes together with $D_1$ as far as the vertical axis in the point $R=1$. The $b_s$ can be earned from stability basins for $\omega$ in the adequate interval for initial value ${\omega }_p=0$. It corresponds to amplitude of $u_{\mathit{Cx}}$ as a testing value for decision about affiliation to $b_{u1}$ or $b_{l1}$ attractiveness.

The interval between $B_1$ and $B_2$ includes the spatial response, see non-trivial amplitude max $\mid u_{\mathit{Cy}}\mid$. The spatial response has a chaotic character, as it has been already outlined in the previous paragraph, when commenting the branch $b_{u2}$.

7. Conclusion

The common physical origin of Lagrangian and Appell-Gibbs approaches has been shown. It originates from the equilibrium of energy-level evolution in time on one side and power supply together with energy dissipation on the other side. Various formulations of this principle lead finally to different variational principles, although they follow from the same minimization of the energy spent to system response portrait. Comparing individual sections of the chapter, we can see that each one of commonly used procedures based on particular energy formulations is preferable for a certain type of problems. It can be concluded that there does not exist a single universal approach which should be recommended.

Some detailed properties of both approaches have been demonstrated in Sections 5 and 6. Both of them discuss non-holonomic problem of the ball movement within the spherical cavity under external excitation. The former one deals with a simple planar problem and shows that the Lagrangian approach is easily applicable to obtain reasonable results as far as a wide parametric discussion, which enable to earn a detailed insight into the system dynamic properties. The latter alternative represents the full space problem with six DOFs and three non-holonomic constraints. Some earlier studies tried to formulate this problem also in Lagrangian style using Lagrange’s multipliers. Finally, it proved that the relevant governing differential system is too complex and does not enable appropriate detailed analysis of dynamic properties of the system. Therefore, the space problem outlined in Section 6 has been formulated using Appell-Gibbs approach. Transparent results have been obtained as needed for practical purposes in a device design and in further study of multi-body system dynamics. Take a note regarding the classification of singular solutions and their applicability for detailed analysis, stability of various regimes of the system under kinematic excitation, transitions among semi-trivial, auto-parametric, chaotic, and other states typical for nonlinear system. Let us add that both 2D and 3D problems have been investigated respecting the full nonlinearity without any simplifications of transcendent functions and thus enabling to study all effects without any limitations in amplitudes.

A certain shortcomings which apparently follow from the knowledge of rotation velocities only (no rotations themselves are calculated) can be disregarded, when displacements have been obtained. The rotation velocities represent mostly satisfactory information. Nevertheless, if rotations are still needed, there exist several variants of a simple differential system (following rotation vector definition) relating velocity and rotation vector components. This system can be subsequently easily solved, when necessary. A hidden complexity of the Lagrangian approach follows from an implicit connection of both parts, which are independent when Appell-Gibbs procedure is applied.

Acknowledgments

The kind support of the Czech Science Foundation project No. 17-26353 J and of the RVO 68378297 institutional support are gratefully acknowledged.