Spatial double physical pendulum with axial excitation: computer simulation and experimental set-up

Physical and mathematical model of a three-dimensional double physical pendulum being coupled by two universal joints is studied. Upper joint of this pendulum is periodically rotated in its axial direction. Damping forces and torques inside joints as well as an influence of the gravitational field and damping forces and torques inside joints are taken into account during derivation of the ODEs of this mechanical system. The work consists of mathematical modelling, computer simulation and parts of the experimental analysis to confirm the numerical simulations. The experimental setup realizes kinetic excitation by an non-constant periodic torque of servomotor, controlled by the computer. Four incremental encoders are mounted on pendulum joints, so the angles of rotation of the pendulum are measured in real-time and analysed by an originally developed acquisition software. Future plans of the mathematical and experimental model are also discussed. Presented exemplary results showed a number of non-linear effects, including chaos, quasi-periodic and periodic dynamics.


Introduction
We consider the mathematical model of a simple threedimensional (3D) double physical pendulums system, under specific conditions. The mechanism is physically excited by a periodic torque in axial direction and linearly damped by each joint. This excitation is applied to the suspension joint of the mechanism. Such an approach (damping and excitation) is often used in many pendulum analysis but rarely the excitation has axial direction.
There are many examples of simple mechanical systems that exhibit complex behavior under certain conditions. Common pendulum, e.g. a clock by Huygens [1], under typical situations, are not enough interesting in term of nonlinear dynamics. It appears here that the conditions are more important than the mechanisms itself. The pendulum became highly fascinating under specific circumstances, e.g. influence of excitation and damping.
The evolution of pendulum analysis started from the measurements and experiments, e.g. famous Foucault's pendulum, 1851 [2,3] or Kater's reversible pendulum [4]. Nowadays pendulum mechanisms are more often used in theoretical or practical mathematical modeling process, for example, to develop more effective vibrations absorption methods [5,6]. There are also advanced research about the pendulum itself, including its experimental identification [7,8] or control algorithms [9].
Multiple pendulum systems are mostly simplified either to planar vibrations [10,11] or they concern only mathematical pendulums [12]. Physical pendulums are examined in their multiple configuration, e.g. [13,14] but also in some simplified forms.
This paper presents the results of numerical computations, as well as the possibilities of applications of the original simulation program. A wide range of nonlinear behavior of the

The Pendulum Model
The system to be considered here is shown in Fig. 1. It is build of two cylindrical-shaped rigid bodies combined by universal joint O 2 and hung on a second universal joint O 1 . This joint is externally driven and it actuates the entire mechanical system axially with either constant or non-constant angular velocity. Influence of the gravitational force is also included. The only damping force that the model contains is inside joints and it is characterized by a simple viscous damping function. Air resistance is neglected due to low velocities of vibrations and relatively large moments of inertia.
Angular positions of each universal joint's shaft have been described by three Euler angles ϕ i , θ i and ψ i , where i corresponds to an index of each joint. The rotation matrices are found, as well as positions of each body centers. Also its linear and angular velocities and energy are defined.
Analytical Wolfram Mathematica computer package has been carried out, during process of derivation of Eq. (1). Full form results are too large and couldn't be simplified enough to show in this paper.
In the study, simple model of viscous damping of joints is assumed in the following form: where M i are corresponding damping torques proportional to the angular velocities. Angular velocity of the axial excitation of first joint follows: where ω 0 is a constant part of velocity (rad/s), q is the amplitude (N m) and Ω states for frequency (rad/s). Unlike experimental setup, the mathematical model does not take into account impacts caused by the mechanical limits of real double spatial pendulum rotation yet. That is why the full comparison of the theoretical model and experimental measurements can be performed for low amplitudes only, smaller than about 90 • . However, this does not prevent to achieve some really interesting simulations results based on the experimental setup parameters (not measurement), without angle limits. Full experimental analysis will be performed after including impacts model to the equations of the motion, increase damping coefficients in the experimental joints to prevent such impact and/or preparing the setup to hold out frequent impacts without the risk of failure.

Numerical computations
Results presented in this paper concern the following fixed parameters (see Fig. 1) presented in Table 1. Values on the left are theoretical ones, chosen for the purpose of exemplary simulations. Two columns on the right presents simulation parameters that corresponds to the measured values of built experimental setup.
According to the User manual of Wolfram Mathematica package, its ODEs solving algorithm is based on higher order Runge-Kutta methods with automatic step control. Results, as well as the plots, are automatically interpolated to any chosen time steps.
Analytically determined ODEs set (1) was simplified and solved numerically for parameters values showed in Table  1. Every first 500 time steps of all numerical computations  were ignored as transient motion and next 500 or more if needed were qualified as significant for the analysis.

Results and analysis
To find globally how much these two systems (example and experiment, see Table 1) are different, 2D maps of maximum numbers of full 360 • rotations were computed (Fig. 2) for the same range of control parameters q and Ω, while ω 0 equals 0 rad/s (pure periodic excitation).
One can see some similarities between them the maps, especially in the shape of the regions where the number of full rotations are low. These maps need a very long time of computing (here about 31 h) but it reveals some general information about the system itself, faster than the bifurcation diagrams.
Using the classical methods of analysis of dynamic systems, the bifurcation diagrams, Poincaré maps and phase plots were produced numerically. In Fig. 3, three nonlinear phenomena are presented. Amplitude of excitation q was set to 12 N m to cover all three regions of different value of full rotations (see Fig. 2).
Similar analysis were performed for the parameters that corresponds to built experimental setup. This time, the system also showed a number of interesting non-linear phenomena, reported in Fig. 4. During the simulation process, it appeared that the system is less sensitive to exhibit nonlinear phenomena than the exemplary one. The reason for this is probably a large moments of inertia of the experimental pendulum, and a small damping of motion. Apart from this, it can be seen here also periodic, quasi-periodic and chaotic behaviour of our investigated system.

Experimental rig
In Fig. 5 one can see the photography and construction scheme of the experimental rig. The original control and acquisition software has been developed, using Java programming language. It can perform in real-time, parallel measurements of values of all four angles of rotation simultaneously. The program can control the parameters of the external excitation and at the same time record each pendulums link position.
The orientation of each pendulum link is measured by four precise incremental encoders and the dedicated PC acquisition card. To transmit signals between rotating pendulums equipped with the encoders and mounting frame, without the risk of the wiring damage, special slip ring is used. It consists of motor-like brushes, which provides a continuous electrical Some preliminary measurements have been performed. Figure 6 shows a time series measurement under constant angular velocity of excitation equals 3.7 rad/s. Presented plot of angular positions of each link shows that as the higher frequency vibrations disappears, all four links oscillate simultaneously. Afterwards, because of the influence of the constant axial excitation, the amplitude of vibration increases gradually to reach its constructional limits, and takes place an impact. This impact sce-   Fig. 7. However, the build experimental setup is not prepared for frequent collisions of its links and has to be modified to prevent its damage.
As showed in Fig. 8, the measurement equipment is enough fast to make some draft phase plots. It should be noted that there is no direct information about the angular velocities of the links, and they have to be calculated numerically. Applied numerical differentiation algorithm makes the plot quite discrete but sufficiently suitable for further generation of bifurcation diagrams or other data analysis.

Conclusions
As it was presented before, double physical pendulum with axial excitation, either in exemplary and experimental configuration, exhibits a wide spectrum of nonlinear effects. Periodic, quasi-periodic and chaotic orbits have been found and discussed, among others. Preliminary measurement data plots have been also showed and shortly discussed. Numerically generated full 360 • rotation maps of the pendulum links show that both configurations reveals some similarity, despite of the fact that its geometrical and mechanical parameters differ significantly. Presented figures and the carried out study confirm that performed analytical and numerical calculations can be used to analyse and simulate multiple pendulum systems and other similar multiple rigid body mechanical systems.
The next step of this project is to extend the mathematical model to include impacts. Additionally, the experimental setup has to be upgraded, to make possibility of modifying the damping coefficient of all pendulums links (ball bearings) and to automatically detect all potential impacts occurrence. Rapid responsible impacts detection is necessary to make a long-time measurements guiding to Poincaré maps and bifurcation diagrams generation. That will make possibility to perform full identifi-cation of built system with no other mechanical limitations.