# The effect of discretization on the numerical simulation of the vibrations of the impacting cantilever beam

Barbara Blazepczyk-Okolewska, Krzysztof Czolczynski, Tomasz Kapitaniak

► To cite this version:

Barbara Blazepczyk-Okolewska, Krzysztof Czolczynski, Tomasz Kapitaniak. The effect of discretization on the numerical simulation of the vibrations of the impacting cantilever beam. Communications in Nonlinear Science and Numerical Simulation, 2010, 15, pp.3073 - 3090. 10.1016/j.cnsns.2009.10.003. hal-01403852

HAL Id: hal-01403852

https://hal.science/hal-01403852

Submitted on 28 Nov 2016

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# The effect of discretization on the numerical simulation of the vibrations of the impacting cantilever beam

Barbara Blaziejczyk-Okolewska, Krzysztof Czolczynski, Tomasz Kapitaniak

Division of Dynamics, Technical University of Lodz, Stefanowskiego 1/15, 90-924 Lodz, Poland

We study the dynamics of a cantilever beam with an unnegligible large mass and with a concentrated mass fixed at the end, which impacts on the base during motion. Generally to model such a system, the finite element method with appropriate number of degrees of freedom has to be employed. However, to analyse some selected aspects of its dynamic behaviour, particularly to predict if the motion with impacts will be periodic, lower-dimen-sional substitutive models with one degree or two degrees of freedom can be employed. The way to determine the parameters of such models and their applicability limits are discussed.

## 1. Introduction

Mechanical systems whose elements impact on one another during operation have been extensively investigated by many researchers. The reason for this interest lies in the fact that motions with impacts exist in a wide variety of engineering applications, particularly in mechanisms and machines with clearances or gaps (e.g., see [1] and the references therein). The physical process during impacts is discontinuous and strongly nonlinear, so the vibro-impact systems can exhibit very rich and complicated dynamic behaviour, e.g., a Feigenbaum scenario (e.g., [2,3]), sudden changes in the chaotic attractor and intermittent to chaos (e.g., [4,5]), Devil's attractors (e.g., [3]), as well as different types of grazing bifurcations (e.g., [6,7]).

The intensive development in investigations of nonlinear behaviours comprises more and more complex mechanical systems with impacts and causes that new scientific teams join the group of researchers involved in dealing with these issues. Vibrations of mechanical systems with impacts have been already discussed in our previous studies (e.g., [4,8-13]). In particular, linear oscillators (e.g., [12]), or sets of linear oscillators (e.g., [8,10,11,13]), whose vibrating motion is disturbed by impacts, have been analysed. These are impacts on a fixed base (e.g., [12]), or impacts of two oscillators (dependent, e.g., [10,13] or independent, e.g., [8,11]) against each other. The investigations have been focused on finding such sets of parameters that characterize the system (mass, stiffness coefficients of springs, coefficients of damping, frequency of external forcing, etc.), for which the motion with impacts is periodic. The investigations are important from the viewpoint of potential applications of vibrating systems with impacts, independently of the fact if motion periodicity is a desirable phenomena or if we want to avoid this periodicity, expecting a chaotic or quasi-periodic motion of the system to occur.
Stability of the periodic motion has also been investigated (e.g., [11,12]) to answer the following issues: does the dis-
turbed periodic motion return to its original form or does it transfer, for instance, into a motion without impacts? How large
are admissible impacts? An effect of inaccuracies in the system manufacturing on its motion periodicity is of great interest as
well: will slight, inevitable differences between the actual and assumed (in the physical model) values of parameters char-
acterizing the system result in the fact that the system will be found unexpectedly beyond the region of the periodic motion?
(see e.g., [11,12]).

In the investigations carried out so far [4,8,10–13] by means of the so-called numerical calculations, physical models with a finite number of degrees of freedom, composed of rigid and heavy masses, connected by a massless spring, have been employed. The investigations have consisted in analysis of the system motion simulated by numerical integrations of the equations of motion. Time series, Poincare maps, bifurcation diagrams and maps of basins of attraction of various coexisting attractors have been prepared for the analysis. It is widely known that such investigations are time consuming, particularly due to a slow stabilization process of disturbed motion of the system, whose parameters are close to the boundary of the periodic motion region.

To avoid these problems we have employed also the analytical method, developed by Peterka [14], for the one degree of freedom system (an oscillator colliding with the unmovable base, e.g., [12]). The method has been generalised to a system with two degrees of freedom [11]. In this method it is possible to find quickly whether there exists a periodic stable solution to the equations of motion for the system under analysis, without numerical integration of these equations. Instead, a possibility of existence of an analytical form of the periodic solution with identical impacts is checked.

In this paper, the investigations so far limited to systems with a low number of degrees of freedom (two degrees at the maximum) with massless springs have been extended to continuous elastic systems (that is to say, systems with considerable, unnegligible mass of elastic elements), that additionally have rigid concentrated masses. The investigations of such systems are related to an important problem of their discretization (discretization of their continuous elements, e.g., vibrating beams, is necessary of course). It is well known that an analysis of continuous systems is limited to the analysis of steady states, and systems themselves cannot have a very complex structure. A question then arises: how to develop discrete models of continuous systems with impacts – how many degrees of freedom should they have and how to define the size of the rigid mass which impacts on the base, etc.?

Reduction in the number of degrees of freedom of the discrete system is also an important issue. One of the reasons is obviously the time-consuming numerical integration of the equations of motion; another one – a possibility of application of the Peterka's method to define the position of regions of the periodic motion. The investigations conducted by the authors so far [9] indicate that an extension of regions of the periodic motion identified with the Peterka's method decreases sharply with an increase in the number of degrees of freedom. This follows from the fact that each impact generates free vibrations of the system, which are a superposition of all eigenmodes of vibrations. As subsequent eigenfrequencies are not their multiplicities in principle (the ratios of eigenvalues are not natural numbers), it is difficult to expect periodicity of such vibrations. In practice, however, due to damping of vibrations in the time intervals between subsequent impacts, the periodicity of vibrations related to, for instance, the two lowest eigenfrequencies only can have a considerable meaning. It seems then that an application of the model with a limited number of degrees of freedom could provide the result closer to reality than an application of the model with a significant number of degrees of freedom.

It is not possible to show a solution to the whole problem of discretization and answer all the above-mentioned questions within the scope of one paper. We restrict thus to a presentation of results of the investigations that were conducted to find an answer to the following question: what are the applicability limits of the model with one degree of freedom or two degrees of freedom as a substitutive system for the system composed of a cantilever beam and a concentrated mass fixed at its end?

We consider the vibrations of the cantilever beam shown in Fig. 1. The beam has a length *l*, mass *m*<sub>b</sub>, inertial momentum of the cross section *l* and is made of the material with the elasticity modulus *E*. One end of the beam is restrained and on the other one the concentrated mass *M* is located. The mass *M* can impact on the base while vibrating. If the system is in the static equilibrium position, the bottom part of the mass *M* (fender) is situated at a distance *d* from the base. A harmonic excitation force of the amplitude *e*η<sup>2</sup> proportional to the square power of its circumferential frequency *η* acts on the mass *M* (e.g., the same amplitude-frequency relation takes place in the case of the force generated by the centrifugal force acting on the rotating rotor).

The investigations were conducted in two stages. In the first stage, for a model shown in Fig. 1 the finite element approximation, hereafter referred to as the FEM model, and the approximations base on one degree of freedom or two degrees of freedom discretizations (hereafter referred to as the 1DOF model or the 2DOF model), shown in Fig. 2(a-c), were developed. In the second stage, various dynamic characteristics (for some selected parameters) determined for the FEM model were

<figure><img src="image_2.png" /><figcaption>Fig. 1. Considered cantilever beam.</figcaption></figure>
<figure><img src="image_3.png" /><figcaption>Fig. 2. Equivalent models of the cantilever beam; FEM model (a), 1DOF model (b); 2DOF model (c).</figcaption></figure>

compared to the respective characteristics for substitute systems, that is to say, for the 1DOF model and the 2DOF model. The results of these comparisons allowed for drawing conclusions on a possibility of using substitute models with a low number of degrees of freedom to model and analyse the motion of the heavy elastic beam with an additional concentrated mass impacting of the base.

## 2. Discrete models

### 2.1. Finite elements method approximation

Fig. 2(a) depicts the discretization of our system generated with the finite element method. Some simple preliminary numerical simulations were conducted to state whether an application of four finite elements to model vibrations of the cantilever beam was sufficient: the values of the first three resonance frequencies of such a model (without the concentrated mass) numerically calculated comply with the analytical results with an accuracy up to three significant digits.

The matrix equation of the system motion in the time intervals between impacts can be written in the following form:

$$ [M]\{\dot{x}\} + [T]\{\dot{x}\} + [K]\{x\} = \{F\} \sin(\eta t). \quad (1) $$

The stiffness matrix $[K]$ takes the form:

$$ [K] = \begin{bmatrix} k_{11}^1 & k_{12}^1 & k_{13}^1 & k_{14}^1 & 0 & 0 & 0 & 0 & 0 & 0 \\ k_{21}^1 & k_{22}^1 & k_{23}^1 & k_{24}^1 & 0 & 0 & 0 & 0 & 0 & 0 \\ k_{31}^1 & k_{32}^1 & k_{33}^1 + k_{11}^2 & k_{34}^1 + k_{12}^2 & k_{13}^2 & k_{14}^2 & 0 & 0 & 0 & 0 \\ k_{41}^1 & k_{42}^1 & k_{43}^1 + k_{21}^2 & k_{44}^1 + k_{22}^2 & k_{23}^2 & k_{24}^2 & 0 & 0 & 0 & 0 \\ 0 & 0 & k_{31}^2 & k_{32}^2 & k_{33}^2 + k_{11}^3 & k_{34}^2 + k_{12}^3 & k_{13}^3 & k_{14}^3 & 0 & 0 \\ 0 & 0 & k_{41}^2 & k_{42}^2 & k_{43}^2 + k_{21}^3 & k_{44}^2 + k_{22}^3 & k_{23}^3 & k_{24}^3 & 0 & 0 \\ 0 & 0 & 0 & 0 & k_{31}^3 & k_{32}^3 & k_{33}^3 + k_{11}^4 & k_{34}^3 + k_{12}^4 & k_{13}^4 & k_{14}^4 \\ 0 & 0 & 0 & 0 & k_{41}^3 & k_{42}^3 & k_{43}^3 + k_{21}^4 & k_{44}^3 + k_{22}^4 & k_{23}^4 & k_{24}^4 \\ 0 & 0 & 0 & 0 & 0 & 0 & k_{31}^4 & k_{32}^4 & k_{33}^4 + k_{sx} & k_{43}^4 \\ 0 & 0 & 0 & 0 & 0 & 0 & k_{41}^4 & k_{42}^4 & k_{43}^4 & k_{44}^4 + k_{stp} \end{bmatrix}. \quad (2) $$
Its components are stiffness matrices of individual finite elements:

$$[K^i] = \begin{bmatrix} k_{11}^i & k_{12}^i & k_{13}^i & k_{14}^i \\ k_{21}^i & k_{22}^i & k_{23}^i & k_{24}^i \\ k_{31}^i & k_{32}^i & k_{33}^i & k_{34}^i \\ k_{41}^i & k_{42}^i & k_{43}^i & k_{44}^i \end{bmatrix} = \frac{EI}{l_e^3} \begin{bmatrix} 12 & 6l_e & -12 & 6l_e \\ 6l_e & 4l_e^2 & -6l_e & 2l_e^2 \\ -12 & -6l_e & 12 & -6l_e \\ 6l_e & 2l_e^2 & -6l_e & 4l_e^2 \end{bmatrix} \quad (3)$$

The support stiffness coefficients $k_{sx}$ and $k_{s\phi}$, whose values are high enough to provide an effect of fixing the beam, are additional elements. Moreover, the symbol E denotes the Young modulus, and the symbol l – a moment of inertia of the beam cross-section.

The inertia matrix has the form:

$$[M] = \begin{bmatrix} m_{11}^i + M & m_{12}^i & m_{13}^i & m_{14}^i & 0 & 0 & 0 & 0 & 0 & 0 \\ m_{21}^i & m_{22}^i & m_{23}^i & m_{24}^i & 0 & 0 & 0 & 0 & 0 & 0 \\ m_{31}^i & m_{32}^i & m_{33}^i + m_{11}^i & m_{34}^i + m_{12}^i & m_{13}^i & m_{14}^i & 0 & 0 & 0 & 0 \\ m_{41}^i & m_{42}^i & m_{43}^i + m_{21}^i & m_{44}^i + m_{22}^i & m_{23}^i & m_{24}^i & 0 & 0 & 0 & 0 \\ 0 & 0 & m_{31}^2 & m_{32}^2 & m_{33}^2 + m_{11}^3 & m_{34}^2 + m_{12}^3 & m_{13}^3 & m_{14}^3 & 0 & 0 \\ 0 & 0 & m_{41}^2 & m_{42}^2 & m_{43}^2 + m_{21}^3 & m_{44}^2 + m_{22}^3 & m_{23}^2 & m_{24}^2 & 0 & 0 \\ 0 & 0 & 0 & 0 & m_{31}^3 & m_{32}^3 & m_{33}^3 + m_{11}^4 & m_{34}^3 + m_{12}^4 & m_{13}^4 & m_{14}^4 \\ 0 & 0 & 0 & 0 & m_{41}^3 & m_{42}^3 & m_{43}^3 + m_{21}^4 & m_{44}^3 + m_{22}^4 & m_{23}^4 & m_{24}^4 \\ 0 & 0 & 0 & 0 & 0 & 0 & m_{31}^4 & m_{32}^4 & m_{33}^4 + M_S & m_{43}^4 \\ 0 & 0 & 0 & 0 & 0 & 0 & m_{41}^4 & m_{42}^4 & m_{43}^4 & m_{44}^4 + B_S \end{bmatrix} \quad (4)$$

This matrix consists of inertia matrices of individual finite elements of masses $m_e$ and lengths $l_e$:

$$[M^i] = \begin{bmatrix} m_{11}^i & m_{12}^i & m_{13}^i & m_{14}^i \\ m_{21}^i & m_{22}^i & m_{23}^i & m_{24}^i \\ m_{31}^i & m_{32}^i & m_{33}^i & m_{34}^i \\ m_{41}^i & m_{42}^i & m_{43}^i & m_{44}^i \end{bmatrix} = \frac{m_b}{420} \begin{bmatrix} 156 & 22l_e & 54 & -13l_e \\ 22l_e & 4l_e^2 & 3l_e & -3l_e^2 \\ 54 & 13l_e & 156 & -22l_e \\ -13l_e & -3l_e^2 & -22l_e & 4l_e^2 \end{bmatrix} \quad (5)$$

The concentrated mass M and the disc of a heavy mass $M_s$ and a high moment of inertia $B_s$ that together with the coefficients $k_{sx}$ and $k_{s\phi}$ yield an effect of fixing the beam, are additional elements.

One of the applied damping models of vibrating systems with many degrees of freedom is the model in which the damping matrix is a sum of two matrices (a detailed description of this model is to be found in, e.g., [15]):

$$[T] = \nu[M] + \zeta[K] \quad (6)$$

The first matrix, proportional to the inertia matrix, is a model of external damping (the coefficient of external damping has been denoted by $\nu$), the second matrix being proportional to the stiffness matrix is a model of internal damping (the coefficient of internal damping is referred to as $\zeta$). In the investigations discussed here, damping was restricted to external damping, assuming that $\zeta = 0.0$.

<figure><img src="image_4.png" /><figcaption>Fig. 3. Eigenfrequencies of the FEM model: $l = 1.0$, $EI = 1/3$.</figcaption></figure>
The structure of vectors of displacements and forces is as follows:

$$ \{x\} = \begin{Bmatrix} x_1 \\ \varphi_1 \\ x_2 \\ \varphi_2 \\ x_3 \\ \varphi_3 \\ x_4 \\ \varphi_4 \\ x_5 \\ \varphi_5 \end{Bmatrix}; \quad \{F\} = \begin{Bmatrix} e\eta^2 \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \end{Bmatrix}. \qquad (7) $$

An impact of the mass M on the base was modelled employing the well-known Newton's law:

$$ v_1^+ = -k_r v_1^- \qquad (8) $$

where $v_1^-$ denotes the velocity before impact and $v_1^+$ stands for the velocity after impact; $k_r$ is the coefficient of restitution. Impacts in the below described substitutive systems with one degree of freedom and two degrees of freedom were modelled in the same way.

To integrate the equations of motion, the Runge-Kutta method was employed; the moments of impacts were established using a method of successive approximations. The information on the dynamic behaviour of the system was provided by time series, phase planes, Poincare maps and bifurcation diagrams made during the integration procedure.

### 2.1.1. DOF discretization

One of the proposed substitutive systems is a system composed of a massless cantilever beam with a concentrated mass fixed at its free end (Fig. 2(b)). The equation of motion of the system in the time intervals between impacts is as follows:

$$ m\ddot{x} + \nu m\dot{x} + kx = e\eta^2 \sin(\eta t), \quad k = \frac{3EI}{l^3}. \qquad (9) $$

The mass m has a fender which can impact on the base during the system motion. If the system is in the static equilibrium position, then the bottom part of the body with the mass m (fender) is situated at a distance d from the base.

### 2.1.2. DOF discretization

The next substitutive system, shown in Fig. 2(c), consists of a massless cantilever beam on which two concentrated masses $m_1$ and $m_2$ are fixed. In this case, the body of the mass $m_1$ has a fender with which it can impact on the base during the system motion. Similarly as in the 1DOF model, when the system is in the static equilibrium position, then the bottom part of the body with the mass $m_1$ (fender) is situated at a distance d from the base.

In the time intervals between impacts, the motion of this system is described with the following equations:

$$ \begin{bmatrix} m_1 & 0 \\ 0 & m_2 \end{bmatrix} \begin{Bmatrix} \ddot{x}_1 \\ \ddot{x}_2 \end{Bmatrix} + \vartheta \begin{bmatrix} m_1 & 0 \\ 0 & m_2 \end{bmatrix} \begin{Bmatrix} \dot{x}_1 \\ \dot{x}_2 \end{Bmatrix} + \begin{bmatrix} KB_{11} & KB_{12} \\ KB_{21} & KB_{22} \end{bmatrix} \begin{Bmatrix} x_1 \\ x_2 \end{Bmatrix} = \begin{Bmatrix} e\eta^2 \\ 0 \end{Bmatrix} \sin(\eta t). \qquad (10) $$

<figure><img src="image_5.png" /><figcaption>Fig. 4. Eigenfrequencies of the FEM model and the 2DOF model; $m_2 = m_{2MAX} = 1.0$ (a); $m_2 = m_{2MAX} = 0.8$ (b).</figcaption></figure>
The stiffness matrix takes the form:

$$ \begin{bmatrix} KB_{11} & KB_{12} \\ KB_{21} & KB_{22} \end{bmatrix} = \frac{1}{W} \begin{bmatrix} 12EI\ell_2^3 & -6EI\ell_2^2(3\ell - l_2) \\ -6EI\ell_2^2(3\ell - l_2) & 12EI\ell_3^3 \end{bmatrix}, \quad (11) $$

where

$$ W = 4\ell_1^3\ell_2^3 - \ell_2^4(3\ell - l_2)^2. \quad (12) $$

The excitation of vibrations in the systems with one degree of freedom and two degrees of freedom takes place analogously as in the system modelled with the finite element method.

It should be mentioned here that our aim has been to derive the 1DOF and 2DOF models which exhibit dynamical behaviour similar to the this of the higher-dimensional reference system (1) so discretization contains parameters which are adjusted to fit the natural frequencies obtained from the FEM model. Alternatively the low-dimensional models can be obtained based on energy principles assuming low number of finite elements [17].

<figure><img src="image_6.png" /><figcaption>Fig. 5. Bifurcation diagrams for $\alpha_2/\alpha_1 = 18.0, \nu = 0.24$; FEM model: $M = 0.79, m_b = 0.21$ (a); 1DOF model: $m = 0.84$ (b); 2DOF model: $m_1 = 0.833, m_2 = 0.167, l_2 = 0.338$ (c).</figcaption></figure>
## 2.2. Eigenfrequencies

The first model investigated was the FEM model with the following parameters: beam length $l = 1.0$, beam rigidity $EI = 1/3$, total mass of the beam and the concentrated mass $m_b + M = 1.0$. The numerical experiments started with calculations of its first two eigenfrequencies $\alpha_1$ and $\alpha_2$. Fig. 3 shows how the values of $\alpha_1$ and $\alpha_2$ alter with a change in the beam mass $m_b$ (and, of course, with a change in the concentrated mass $M = 1 - m_b$). For the value of $m_b = 0$ (a massless beam with the concentrated mass $M = 1$), the frequency $\alpha_2 \to \infty$, whereas $\alpha_1 = 1$. For $m_b = 1$ (a heavy beam without a concentrated mass at its end), the eigenfrequencies equal to $\alpha_1 = 2.029$ and $\alpha_2 = 12.731$ were obtained, employing four finite elements. These values agree with the values from the analytical formulas (commonly published, see, e.g., [16]), equal to $\alpha_1 = 2.0297$ and $\alpha_2 = 12.721$ - differences arise only in the fourth significant digit. An application of a lower finite element number yields the values $\alpha_1 = 2.029$ and $\alpha_2 = 12.75$ for three finite elements and $\alpha_1 = 2.034$ and $\alpha_2 = 12.82$ for two finite elements, respectively. During the former investigations [9] it was stated that the ratio of eigenfrequencies $\alpha_2/\alpha_1$ exerted a considerable influence on the system motion periodicity: the ranges of the circumferential frequency of the excitation force for which the system motion is periodic are wider when a value of the ratio $\alpha_2/\alpha_1$ is close to an even number than if a value of the ratio $\alpha_2/\alpha_1$ is close to an odd number. Therefore, the values of $\alpha_2/\alpha_1$ are shown additionally in Fig. 3 as well.

Fig. 4(a and b) presents a comparison of the values of $\alpha_1$ and $\alpha_2$ for the FEM model with the values of $\alpha_1$ and $\alpha_2$ obtained for the 2DOF model with the parameters: beam length $l = 1.0$, beam rigidity $EI = 1/3$, total values of the concentrated masses $m_1 + m_2 = 1.0$. Apart from the value of masses $m_1$ and $m_2$, the parameter $l_2$ (see Fig. 2b) is a decisive parameter as far as values of eigenfrequencies are concerned. Fig. 4(a) shows the values of $\alpha_1$ and $\alpha_2$ for the 2DOF model (thin lines) calculated for four different values of this parameter: $l_2 = 0.2$, $l_2 = 0.3$, $l_2 = 0.4$, and $l_2 = 0.5$, correspondingly. In Fig. 4(a) on the horizontal axis, the value of the mass $m_2$ alters from $m_2 = 0.0$ up to $m_2 = 1.0$. The values of eigenfrequencies for the FEM model (thick lines) are presented as well. It can be easily seen that in the wide range $0 < m_2 < 0.5$, the values $\alpha_1$ of the 2DOF model and the FEM model are almost identical; at $l_2 = 0.3$, the values $\alpha_2$ of the 2DOF model are close to $\alpha_2$ for the FEM model. Considerable differences between the values of both $\alpha_1$ and $\alpha_2$ can be seen in the right-hand side of the figure. It follows from the fact that $m_2 \to 1.0$ means that $m_1 \to 0.0$, and thus a return to the 1DOF model, that is to say, $\alpha_2 \to \infty$, whereas for the FEM model $m_b \to 1.0$ means that $\alpha_2$ tends to the value $\alpha_2 = 12.731$. Hence, the maximum value of the mass $m_2$ has to be chosen. This value is denoted by $m_{2MAX}$, whereas the maximum value on the horizontal axis $m_2 = m_{2MAX}$ will correspond to the maximum value $m_b = 1.0$ for the FEM model. During further experiments, it was decided that the optimal value was $m_{2MAX} = 0.8$.

Fig. 4(b) shows once again the values $\alpha_1$ and $\alpha_2$ for the FEM model in the range from $m_b = 0$ up to $m_b = 1$ and the values $\alpha_1$ and $\alpha_2$ for the 2DOF model in the range from $m_2 = 0$ to $m_{2MAX} = 0.8$; the parameter defining the position of the mass $m_2$ was

<figure><img src="image_7.png" /><figcaption>Fig. 6. Bifurcation diagrams for $\alpha_2/\alpha_1 = 18.0, \nu = 0.08$; FEM model: $M = 0.79, m_b = 0.21$ (a); 2DOF model: $m_1 = 0.833, m_2 = 0.167, l_2 = 0.338$ (b).</figcaption></figure>
given two close values: $l_2 = 0.338$ and $l_2 = 0.315$. As results from the analysis of this figure, a slight change in $l_2$ allows one to obtain the values $\alpha_1$ and $\alpha_2$ for the 2DOF that comply with the FEM model.

### 3. Comparison of the dynamic behaviour of the FEM model and the substitutive 1DOF and 2DOF models

During further numerical investigations, the dynamic behaviour presented on bifurcation diagrams prepared for the FEM models with selected values of the beam mass was compared with the dynamic behaviours obtained on bifurcation diagrams made for the 1DOF model and the 2DOF model. The value of the mass $m$ for the 1DOF model was calculated so that its eigen-frequency $\alpha$ was equal to the fundamental eigenfrequency $\alpha_1$ of the FEM model. As $3EI = 1.0$ and $l = 1.0$, a simple relationship is obtained:

$$ \alpha_1^2 = \alpha^2 = \frac{3EI}{ml^3} \Rightarrow m = \frac{1}{\alpha_1^2} \qquad (13) $$

The parameters of the 2DOF model were selected in such a way that its both eigenfrequencies were equal to the corresponding eigenfrequencies of the FEM model, and the total mass $m_1 + m_2$ was equal to the mass $m_b + M$. To fulfil these conditions, suitable values of the parameters $m_2$ and $l_2$ that defined the size and position of the mass $m_2$ were assumed – see Fig. 4(b).

The investigations started with a system in which 79% of the mass was the concentrated mass at the beam end. It was forecast that in such a case it would be easier to select substitutive systems with similar dynamic behaviours (in particular, the characteristics of the bifurcation diagram that allows for investigations of the system dynamics in a wide range of the control parameter). Fig. 4(b) provides the information that for the FEM model with the concentrated mass at the beam end $M = 0.79$ and the beam mass $m_b = 0.21$, the eigenfrequencies are equal to $\alpha_1 = 1.092$ and $\alpha_2 = 19.72$, and the ratio of frequencies $\alpha_2/\alpha_1 = 18.0$. The damping coefficient $\nu = 0.24$ (which corresponds to the value of the logarithmic decrement $\Delta = \ln(2)$ for the system with 1DOF) has been assumed. Such a high value of the damping coefficient was chosen to minimize an effect of vibrations with the frequency $\alpha_2$ on the periodicity of motion with impacts.

Fig. 5(a) shows a bifurcation diagram made for the above-described FEM model. The value $e = 1.0$ was assumed, i.e., vibrations were excited with the harmonic force $F(t) = \eta^2 \sin(\eta t)$, whose circumferential frequency varied from $\eta = 3.0$ to

<figure><img src="image_8.png" /><figcaption>Fig. 7. Comparison of the methods for the 2DOF model for $m_1 = 0.833$, $m_2 = 0.167$, $\alpha_2/\alpha_1 = 18.0$, $\nu = 0.0$; regions of stable solutions of the Peterka's map (a); bifurcation diagram obtained by means of an integration of motion equations with the Runge-Kutta method: $l_2 = 0.338$ (b).</figcaption></figure>
$\eta = 10.0$. The graduation on the vertical axis of the plot represents a value of displacement of the concentrated mass $x_1$ at the time instants defined by the equation $\sin(\eta\tau) = 1.0$. Impacts on the fixed base took place during motion; the coefficient of restitution $k_r = 0.6$ and the distance $d = 0.0$ were assumed. The diagram shows how the character of the system motion depends on the frequency value of the excitation force $\eta$. Three distinct regions of periodic vibrations can be observed: with period 2 ($3.6 < \eta < 5.1$), with period 3 ($5.8 < \eta < 7.2$), and with period 4 ($8.0 < \eta < 9.5$). For the remaining values of $\eta$, we can observe nonperiodic vibrations: *quasi-periodic* or *chaotic*.

The displacements $x$ of the substitutive mass of the value $m = 0.839$ and the frequency $\alpha = 1.092$ for the 1DOF model have been presented on the vertical axis of the bifurcation diagram denoted as Fig. 5(b). The value of the damping coefficient $\nu$, the way of excitation of vibrations and the parameters $k_r$ and $d$ were assumed identical as for the FEM model.

It follows from Fig. 4(b) that the substitutive 2DOF model has eigenfrequencies $\alpha_1 = 1.092$ and $\alpha_2 = 19.72$ if $m_1 = 0.833, m_2 = 0.167$ and $l_2 = 0.338$. Fig. 5(c) presents a bifurcation diagram made for this model, when the way of excitation of vibrations and the parameters $\nu, k_r$ and $d$ are assumed identical as for the FEM model. The graduation on the vertical axis of the diagram denotes displacements $x_1$ of the concentrated masses $m_1$.

A comparison of the bifurcation diagrams for the FEM model and the respective 1DOF model and the 2DOF model shows a very good compatibility of both the ranges of periodic and nonperiodic motion, as well as amplitudes of vibrations – which,

<figure><img src="image_9.png" /><figcaption>**Fig. 8.** Bifurcation diagrams for $\alpha_2/\alpha_1 = 8.0, \nu = 0.33$; FEM model: $M = 0.293, m_b = 0.707$ (a); 1DOF model: $m = 0.461$ (b); 2DOF model: $m_1 = 0.452, m_2 = 0.548, l_2 = 0.315$ (c).</figcaption></figure>
as has been mentioned, has been expected for models with a heavy concentrated mass at the beam end. On the bifurcation diagrams (Fig. 5(a) and (c)), the values of displacement $x_1$ decrease in the vicinity of the excitation frequency $\eta = 4.37$, $\eta = 6.55$ and $\eta = 8.74$. This phenomenon is caused by an increase in the amplitude of free vibrations of the frequency $\alpha_2$, generated by impacts. A more detailed explanation of this issue is to be found later in this paper.

A decrease in the value of the damping coefficient from $v = 0.24$ to $v = 0.08$ for vibrations of the frequency $\alpha_1$ is not followed by differences arising between the bifurcations diagrams made for the FEM model (Fig. 6(a)) and the 2DOF model (Fig. 6(b)): one cannot see any differences in the position of ranges of the periodic motion and the amplitudes of periodic vibrations.

Fig. 7(a) depicts a plot drawn by means of the Peterka's method (a detailed description of the method is to be found, for instance, in [9] or [14]) of the undamped ($v = 0.0$) 2DOF model with the masses $m_1 = 0.833$ and $m_2 = 0.167$. The frequency $\eta$ and the length $l_2$ are parameters on the horizontal and vertical axis, respectively. One can see the regions of $\eta$ and $l_2$, marked in the gray scale, in which the motion of the system under analysis is periodic with period 2, 3 and 4, on this plot. Moreover, a dependence of the motion periodicity, described in [9], on the value of the ratio $\alpha_2/\alpha_1$, which in turn depends on the value of $l_2$, can be seen as well. Fig. 7(b) shows displacements of the mass $m_1$ of the 2DOF model drawn on the basis of the numerical integration of the equations of motion with the Runge-Kutta method for the selected length $l_2 = 0.338$. A comparison of both the plots proves that the ranges of frequency $\eta$ in which the motion is periodic, identified with both the methods, are the same. This fact confirms the applicability and meaning of the Peterka's method as a method which allows one to find quickly a dependence between the motion parameters of the system with impacts and its motion periodicity. A comparison of Fig. 5(c) ($v = 0.24$), Fig. 6(b) ($v = 0.08$) and Fig. 7(b) ($v = 0.0$) shows that damping increases the ranges of frequency $\eta$, in which the motion is periodic, and thus the Peterka's method is more restrictive and – hence – safer. At low damping, a decrease in the displacement $x_1$ in the vicinity of $\eta = 6.55$ and $\eta = 8.74$ (due to the domination of vibrations with the frequency $\alpha_2$) results in additional impacts on the base and a nonperiodic motion – one can see blurred columns of points instead of two or three points on the bifurcation diagram (see Fig. 7(b)).

## 4. Influence of the beam mass

The subsequent experiments were devoted to systems with a heavier and heavier beam mass. An increase in the beam mass is followed by a decrease in the value of ratio $\alpha_2/\alpha_1$. During the analysis of bifurcation diagrams with a decreasing value of $\alpha_2/\alpha_1$, significant differences (that is to say, differences in the values of amplitudes) between the diagram for the FEM

<figure><img src="image_10.png" /><figcaption>Fig. 9. Bifurcation diagrams for $\alpha_2/\alpha_1 = 8.0, v = 0.11$; FEM model: $M = 0.293, m_b = 0.707$ (a); 2DOF model: $m_1 = 0.452, m_2 = 0.548, l_2 = 0.315$ (b).</figcaption></figure>
model and the corresponding diagrams for the substitutive systems were noticed no sooner than for the system of $\alpha_2/\alpha_1 = 8.0$ and $v = 0.33$ (the value $v = 0.33$ corresponds to the logarithmic decrement of damping $\Delta = \ln(2)$), with the following parameters: $M = 0.293$, $m_b = 0.707$ for the FEM model (see Fig. 8(a)), $m = 0.461$ for the 1DOF model (see Fig. 8(b)), and $m_1 = 0.548$, $m_2 = 0.452$, $l_2 = 0.318$ for the 2DOF model (see Fig. 8(c)). It is worth noticing that a comparison of the bifurcation diagrams shown in Fig. 8(a-c), despite the above-described amplitude divergence, still shows a good coexistence of both the ranges of periodic and nonperiodic motion. It should be added that there is also conformity of values of vibration amplitudes for the 1DOF model (Fig. 8(b)) and the 2DOF model (Fig. 8(c)). It was checked that there were also such conformities for the frequency $\eta$ values higher than the value of the second eigenmode. In the vicinity of the frequency $\eta \approx 5.9$, for

<figure><img src="image_11.png" /><figcaption>Fig. 10. Time series of the FEM model for $M = 0.293$, $m_b = 0.707$, $\alpha_2/\alpha_1 = 8.0$, $v = 0.11$; $\eta = 5.6$ (a); $\eta = 5.86$ (b); $\eta = 5.93$ (c); $\eta = 5.53$ (d).</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Fig. 11. Details of the time series of the FEM model for $M = 0.293$, $m_b = 0.707$, $\alpha_2/\alpha_1 = 8.0$, $v = 0.11$, $\eta = 5.86$.</figcaption></figure>
which the time interval between impacts equal to a doubled period of the forced vibrations $2T_η$ is equal to the free vibration period $T_{α2}$, distinct differences in the value of amplitude and phase (a phase shifting) were observed between the FEM model and the 2DOF model (the same behaviours but for different values of the system parameters and the value $η$ can also be observed in Fig. 5(a-c)). Here we can see a transition through superharmonic resonances. A similar phenomenon can be observed at $η = 8.83$, for which the time interval between impacts is equal to $3T_η = 4T_{α2}$. On the other hand, a local decrease of $x_1$ in the neighbourhood of the frequency $η = 5.53$ is caused by free vibrations with the frequency $α_3$. It was calculated that for the investigated FEM model of $α_2/α_1 = 8.0$, the frequency $α_3 = 35.92$. It means that for the $η = 5.53$, the equality $2T_η = 13T_{α3}$ holds. A different behaviour of both the systems is caused by various mass distributions: the concentrated mass $M = 0.293$ is lighter than the mass $m_1 = 0.452$; the beam mass $m_b = 0.707$ is heavier than the mass $m_2 = 0.548$. Moreover, the mass $m_2$ is fixed at one third of the beam length, whereas the beam mass centre in the FEM model lies (obviously) at its half-length. A detailed analysis of the way these phenomena originate is presented in the further part of this paper for systems with lower damping.

**5. Influence of damping**

As has been mentioned above, slight damping can limit the applicability of models with a low number of degrees of freedom in predicting the periodicity of motion of continuous systems: for the FEM model with the ratio $α_2/α_1 = 8.0$ and the frequency $η = 6.64$, a system chaotic motion appears as a result of a decrease in damping from $v = 0.33$ to $v = 0.11$ (damping $v = 0.11$ corresponds to $Δ = \ln(1.25)$). Fig. 9(a) shows a bifurcation diagram made for this model of the following parameters: $M = 0.293$, $m_b = 0.707$, $α_1 = 1.473$, $α_2 = 11.77$, $v = 0.11$. A decrease in damping has caused that only narrow windows of a periodic motion can be seen in Fig. 9(a), if compared to Fig. 8(a). Fig. 9(b) presents a bifurcation diagram for the respective 2DOF model of the parameters as follows: $m_1 = 0.452$, $m_2 = 0.548$, $l_2 = 0.315$, $v = 0.11$ and, as one can see, a decrease in damping results only in a slight narrowing of ranges of a periodic motion. Moreover, in the vicinity of $η = 5.9$ and $η = 8.83$, additional impacts (see Fig. 8(c)) that are followed by a loss of motion periodicity can be observed for the FEM model.

Fig. 10 depicts time series for the FEM model of the parameters: $M = 0.293$, $m_b = 0.707$, $α_2/α_1 = 8.0$ and $v = 0.11$ (the parameters as on the bifurcation diagram in Fig. 9(a)), for four values of frequency: $η = 5.6, 5.86, 5.93$ and $5.53$. The displacement $x_1$ of the mass $M$ (denoted by 1) and the displacements $x_2, x_3$ and $x_4$ of the nodes connecting the finite elements

<figure><img src="image_13.png" /><figcaption>**Fig. 12.** Comparison of the methods for the 2DOF model for $m_1 = 0.452$, $m_2 = 0.548$, $α_2/α_1 = 8.0$, $v = 0.0$; bifurcation diagram obtained by means of an integration of motion equations with the Runge-Kutta method: $l_2 = 0.315$ (a); regions of stable solutions of the Peterka's map (b).</figcaption></figure>
versus time (denoted by **2, 3, 4**) can be seen in the figure. The time is represented by the number *N* of periods of external forcing, related to time by the relationship *N* = *ηt*/(2π. It is easily seen that at the frequency *η* = 5.6 (Fig. 10(a)), the first eigenmode of vibrations and forced vibrations are the main components of the beam deflection line (and hence, of the displacements of nodes depicted in the figure). A contribution of the second component – the second eigenmode of vibrations – is almost invisible in the displacement *x*₁, however, it can be easily observed in the displacement *x*₄ of the node that is situated closest to the fixed beam end. The system performs the period 2 motion. Fig. 10(b) shows time series of vibrations for *η* = 5.86, at which 2*T*<sub>*η*</sub> = 4*T*<sub>2x2</sub>. For this value of *η*, a contribution of the second eigenmode of vibrations is predominant. The details of the time series of displacements have been enlarged in Fig. 11. Just before impact, the whole beam displaces with the concentrated mass towards the base (downwards), and moreover values of displacements of nodes are lower than the displacement of the concentrated mass, which means that the beam is deflected towards the base. It has also been stated (which can be seen in the figure as well) that the velocity of the concentrated mass before impact is lower than at *η* = 5.6. After bouncing from the base, the concentrated mass displaces upwards, whereas the beam, which has not impacted on anything, displaces still downwards. If the beam is heavy enough in comparison to the concentrated mass (and this is the case here: *m*<sub>b</sub> = 0.707, *M* = 0.293), then the beam draws the concentrated mass behind itself, causing thus an additional impact. And this is the reason why the system motion becomes nonperiodic. Just after the resonance, at *η* = 5.93 (Fig. 10(c)), the beam is deflected upwards at the instant of impact – displacements of its nodes are higher than the displacement of the concentrated mass – besides, node velocities are low. Such a configuration of the beam makes it possible for the concentrated mass to separate from the base after impact – additional impacts do not take place and the system performs the period 2 motion again. Summing up, in Fig. 10 we can see reasons of the periodicity loss of motion of the FEM model: additional impacts appear when the second eigenmode of vibrations is predominant in the beam deflection line and the beam is heavy enough if compared to the concentrated mass. In Fig. 9(b), showing a bifurcation diagram of the 2DOF model, we do not observe this periodicity loss due to a different mass distribution, which has already been mentioned in comments on Fig. 8. It has also been mentioned there that the stability loss of motion of the FEM model may be caused by free vibrations of the subsequent eigenfrequency *α*<sub>3</sub>, generated during impacts: at *η* = 5.53, an equality 2*T*<sub>*η*</sub> = 13*T*<sub>2x3</sub> holds. Fig. 10(d) depicts time series of vibrations for this value of the excitation frequency. One can see that the third eigenmode of vibrations is predominant and that it results in additional impacts – the system motion is nonperiodic. An increase in the excitation frequency up to *η* = 5.6 (Fig. 10(a)) causes that the contribution of the third mode decreases, which makes the motion periodic again.

<figure><img src="image_14.png" /><figcaption>**Fig. 13.** Time series of the 2DOF model for *m*₁ = 0.452, *m*₂ = 0.548, *l*₂ = 0.315, *α*₂/*α*₁ = 8.0, *v* = 0.0; *η* = 5.86 (a); *η* = 5.871 (b); *η* = 5.874 (c); *η* = 5.9 (d).</figcaption></figure>
A question arises then: can additional impacts be observed in the 2DOF model as well? In looking for an answer, the damping in the 2DOF model was decreased to zero ($v = 0.0$) and Fig. 12 was drawn. On the bifurcation diagram in Fig. 12(a), the regions of periodic motion that can be compared to the regions of periodic motion obtained with the Peterka's method (Fig. 12(b)) can be seen. We observe here a very good agreement of the results again. The most interesting, however, is the fact that for the first time we can see the motion periodicity loss in the range $5.86 < \eta < 5.875$, that is to say, in the vicinity of the resonance $2T_{\eta} = 4T_{x2}$, where four periods of free vibrations with the second eigenmode occur in the time interval between impacts. This phenomenon can be observed in the FEM model with a different mass distribution even at strong damping. Fig. 13 shows a time series of vibrations for the 2DOF model for four values of frequencies: $\eta = 5.86, 5.871, 5.874$ and $5.9$. First of all, let us draw attention to the fact that before the resonance ($\eta = 5.86$, Fig. 13(a)), the vibrations $x_1$ and $x_2$ of the masses $m_1$ and $m_2$ are in phase, whereas after the resonance ($\eta = 5.9$, Fig. 13(d)) are in anti-phase. At the frequency $\eta = 5.86$, the mass $m_1$ velocity after impact is so high that the mass $m_2$, still displacing towards the base, cannot draw $m_1$ behind itself towards the base – as a result, an additional impact does not occur. An increase in the excitation frequency changes the phase shifting between the displacements of both masses. Consequently, the mass $m_1$ velocity decreases after impact and the mass $m_2$ is able to change the direction of the mass $m_1$ motion, leading to an additional impact ($\eta = 5.871$, Fig. 13(b)), and – with a further increase in the excitation frequency – even to a series of impacts

<figure><img src="image_15.png" /></figure>

<figure><img src="image_16.png" /></figure>

<figure><img src="image_17.png" /><figcaption>**Fig. 14.** Bifurcation diagrams for $\alpha_2/\alpha_1 = 7.0$ and $\Delta = \ln(2)$ for all models; FEM model: $M = 0.171, m_b = 0.829, v = 0.36$ (a); 1DOF model: $m = 0.37, v = 0.275$ (b); 2DOF model: $m_1 = 0.36, m_2 = 0.64, l_2 = 0.298, v = 0.275$ (c).</figcaption></figure>
occurring one after another, which in real systems can be interpreted an instantaneous continuous contact of the mass $m_1$ with the base (Fig. 13(c)). A further change in the phase shifting restores the period 2 motion of the system.

## 6. System with a very heavy beam

Next, the FEM model of the following parameters: $M = 0.171, m_b = 0.829, \alpha_1 = 1.646, \alpha_2 = 11.52, \alpha_3 = 34.12, \alpha_2/\alpha_1 = 7.0, v = 0.36$, whose bifurcation diagram is presented in Fig. 14(a), and the corresponding 1DOF model ($m = 0.369, v = 0.275$) and the 2DOF model ($m_1 = 0.36, m_2 = 0.64, l_2 = 0.298, v = 0.275$), whose bifurcation diagrams are in Figs. 14(b) and 14(c), respectively, were investigated. The present approach to the selection of parameters makes it possible for various 1DOF, 2DOF and FEM models to assume the parameter $v$ in such a way that the velocity with which free vibrations initiated by the boundary conditions proportional to the first eigenmode of vibrations vanish is the same (and thus, the logarithmic decrement of damping $\Delta = \ln(2)$ for the appropriate set of models).

First of all, let us focus our attention on a considerable decrease in the regions of periodic motion in comparison to the respective ranges marked in Figs. 5(a) and 8(a), drawn for the FEM models with the frequency ratio equal to $\alpha_2/\alpha_1 = 18$ and $\alpha_2/\alpha_1 = 8$, and with the damping coefficient value $v$ corresponding to the logarithmic decrement of damping $\Delta = \ln(2)$ for vibrations with the frequency $\alpha_1(v = 0.24)$ in Fig. 5(a) and $v = 0.33$ in Fig. 8(a)). We guess that the reason are again free vibrations of the frequencies $\alpha_2$ and $\alpha_3$, that will lead to additional impacts in a wider and wider range of the frequency $\eta$ at a decreasing concentrated mass $M$. Fig. 15 shows a time series of vibrations for the FEM model for four values of frequency $\eta = 5.93, 5.76, 5.686$ and $5.25$. At $\eta = 5.93$ (Fig. 15a), we can see the period 2 motion with one impact per period. Fig. 14(a) indicates that the region $\eta$ of periodic motion must be very narrow. A decrease in the excitation frequency up to $\eta = 5.76$ (Fig. 15(b)), at which $2T_\eta = T_{\alpha2}$, brings about an additional impact, whereas a further decrease in the excitation frequency up to $\eta = 5.686$ (Fig. 15(c)), for which $2T_\eta = 12T_{\alpha3}$, results in two additional impacts. Similarly, after decreasing the excitation frequency to $\eta = 5.25$ (Fig. 15(d)), for which $2T_\eta = 13T_{\alpha3}$, three additional impacts appear. As a consequence of these additional impacts, the system motion is nonperiodic, that looks like the period 2 motion.

The next issue is a compatibility of the regions of periodic motion that are observed on the bifurcation diagrams presented in Fig. 14(b) and (c) with the results obtained with the Peterka's method. Fig. 16(a) shows the regions of periodic

<figure><img src="image_18.png" /><figcaption>**Fig. 15.** Time series of the FEM model for $M = 0.171, m_b = 0.829, \alpha_2/\alpha_1 = 8.0, v = 0.36; \eta = 5.93$ (a); $\eta = 5.76$ (b); $\eta = 5.686$ (c); $\eta = 5.25$ (d).</figcaption></figure>
<figure><img src="image_19.png" /><figcaption>**Fig. 16.** Comparison of the methods for the 2DOF model for $m_1 = 0.360$, $m_2 = 0.640$, $\alpha_2/\alpha_1 = 7.0$; regions of the stable solutions of the Peterka's map with $v = 0.0$ (a); bifurcation diagrams obtained by means of an integration of motion equations with the Runge-Kutta method: $l_2 = 0.298$ and $v = 0.01$ (b) and $v = 0.0$ (c).</figcaption></figure>

motion determined with this method for the above-described 2DOF model, as a function of the frequency $\eta$ and the length $l_2$. As has been mentioned in [9], the value of ratio $\alpha_2/\alpha_1 = 7$, i.e., it is an odd number and it denotes very narrow regions of periodic motion obtained with the Peterka's method, which can be seen in Fig. 16(a) for $l_2 = 0.298$. They are much narrower than the regions marked in Fig. 14(c), plotted on the assumption of strong damping $v = 0.275$. When damping is decreased to a low value of $v = 0.01$, on the bifurcation diagram (Fig. 16(b)) one can see new regions of nonperiodic motion caused by additional impacts, but still the regions of periodic motion are much wider than it follows from the calculations conducted with the Peterka's method. A conclusion can be drawn that even slight damping (that always exists in real systems) decreases the amplitude of free vibrations with the frequency $\alpha_2$ to such an extent that these vibrations do not give rise to additional impacts of the mass $m_1$ on the base and the effect of decreasing the regions of periodic motion caused by the ratio of eigenfrequencies close to an odd number is not observed. When damping is decreased to $v = 0.0$ (Fig. 16(c)), a compatibility
of positions of regions of *the period 1 motion* (in the vicinity of $\eta = 4$) and of *the period 2 motion* (in the vicinity of $\eta = 7.5$) is to be observed. No regions of periodic motion in the vicinity of $\eta = 3.3, \eta = 5.8, \eta = 6.7$ and $\eta = 9.0$ can be seen on the bifurcation diagram. The reason can lie in a high sensitivity of the motion to changes in the system parameters, e.g., for $\eta = 6.7$, a change in $l_2$ from $l_2 = 0.298$ to $l_2 = 0.295$ causes that the calculations with the Peterka's method indicate the nonperiodicity of motion.

## 7. Conclusions

Substitutive systems with one degree of freedom, and especially with two degrees of freedom, can be employed in the investigations of periodicity of motion with impacts of simple vibrating systems with a substantial mass of elastic elements. We owe this possibility to the fact that forced vibrations and the first eigenmode of these vibrations are the dominant component of the system geometrical configuration during the motion with impacts.

The higher the similarity (as regards amplitudes and positions of regions of periodic motions) of bifurcation diagrams of the model based on the finite element method (FEM model) and simple models with one degree of freedom (1DOF model) and with two degrees of freedom (2DOF model), the lower the mass of elastic elements. Then, all models maintain a similar mass distribution. In the case of systems with heavy elastic elements, free vibrations of higher frequencies (an effect of vibrations with the frequencies $\alpha_2$ and $\alpha_3$ has been investigated) result in additional impacts and, consequently, in a loss of motion periodicity. In substitutive systems with massless elastic elements, free vibrations with higher frequencies have a lower contribution in the deflection line (if they occur at all due to a limited number of degrees of freedom) because of a different mass distribution. Hence, the regions of periodic motion can be wider in these systems. The less considerable the differences are, the stronger the damping of vibrations is. Thus, the heavier the beam mass, the stronger the damping has to be to make the motion of the FEM model and the corresponding 1DOF and 2DOF models periodic in the same ranges of the frequency $\eta$. On the other hand, from the optimistic point of view, the above-described irregularities and differences occurred only for the system in which the beam mass was equal to 70% of the total system mass.

It has been found that the basic condition the substitutive 1DOF model has to satisfy is maintaining the same fundamental eigenfrequency as the model based on the finite element method. The substitutive 2DOF model must have eigenfrequencies equal to the first and second eigenmode of the FEM model, and besides it must keep its total mass (that is to say, $m_1 + m_2 = m_b + M$). Numerical experiments with various substitutive 2DOF models, prepared according to other criteria, have been conducted. For instance, the FEM model and the substitutive 2DOF model have to have the same values of the eigenvalues $\alpha_1$ and $\alpha_2$ and the mass centre position. Let us notice that this newly formulated condition can replace only the currently used one (and not to be parallel to it, i.e., the total mass of the 2DOF model can be, or even must be different from the mass of the FEM model). These experiments have provided bifurcation diagrams that differ evidently from the diagrams for FEM models and have been concluded with the statement that the 2DOF model maintaining eigenfrequencies and mass of the FEM model is the best (that is to say, the authors cannot develop any better one).

It was stated during the numerical experiments that if a periodic motion is desirable from the viewpoint of the application of vibrating systems, then only a motion with one impact per period can be employed in the engineering practice. An appearance of additional impacts results in a nonperiodic motion. It was also found that a periodic motion with two or three impacts (that occur one after another) per one period of the system motion can occur in some ranges of the frequency of the external excitation $\eta$. However, these intervals are so narrow (of the magnitude of 0.0001) that they can be neglected from the practical point of view.

In theory, there is of course a possibility to generate additional impacts by free vibrations of even higher frequencies: $\alpha_4, \alpha_5$, and so on. However, due to the self-evident reason (width of the range $\eta$), the search for such instances of the motion periodicity loss has been found purposeless.

## Acknowledgment

This study has been supported by the Polish Department for Scientific Research (DBN) under Project No. N501 037 31/2508.

## References

[1] Jerrelind J, Stensson A. Nonlinear dynamics of parts in engineering systems. Chaos Solitons Fractals 2000;11:2413-28.

[2] Isomaki H, Boehm J, Raty R. Devil's attractors and chaos of a driven impact oscillator. Phys Lett A 1985;107:343-6.

[3] Thompson JMT, Ghaffari R. Chaos after period-doubling bifurcations in the resonance of an impact oscillator. Phys Lett A 1982;91:296-9.

[4] Blazejczyk-Okolewska B, Czolczynski K, Kapitaniak T, Wojewoda J. Chaotic mechanics in systems with impacts and friction. World scientific series in nonlinear science series A. Singapore: World Scientific Publishing Co.; 1999.

[5] Peterka F, Vacik J. Transition to chaotic motion in mechanical systems with impacts. J Sound Vibration 1992;154:95-115.

[6] Chin W, Ott E, Nusse HE, Grebogi C. Grazing bifurcations in impact oscillators. Phys Rev E 1994;50:4427-44.

[7] Nordmark AB. Non-periodic motion caused by grazing incidence in an impact oscillator. J Sound Vibration 1991;145:279-97.

[8] Blazejczyk-Okolewska B, Brindley J, Czolczynski K, Kapitaniak T. Antiphase synchronization of chaos by noncontinuous coupling: two impacting oscillators. Chaos Solitons Fractals 2001;12:1823-6.

[9] Blazejczyk-Okolewska B, Czolczynski K, Kapitaniak T. Dynamics of two degree of freedom cantilever beam with impacts. Chaos Solitons Fractals 2007. doi:10.1016/j.chaos.2007.09.097.
[10] Czołczyński K, Kapitaniak T. Influence of the mass and stiffness ratio on a periodic motion of two impacting oscillators. Chaos Solitons Fractals 2003;17:1–10.

[11] Czołczyński K, Kapitaniak T. On the existence of a stable periodic solution of two impacting oscillators with damping. Int J Bifurcation Chaos 2004;14:3931–47.

[12] Czołczyński K, Kapitaniak T. On the existence of a stable periodic solution of an impacting oscillator with two fenders. Int J Bifurcation Chaos 2004;14:3115–34.

[13] Peterka F, Blazejczyk-Okolewska B. Some aspects of the dynamical behavior of the impact damper. J Vibration Control 2003;11:459–79.

[14] Peterka F. An investigation of the motion of impact dampers, paper I, II, III. Strojnicku Casopis XXI, c.5; 1971.

[15] Bathe Jr K, Wilson EL. Numerical methods in finite element analysis. Englewood Cliffs: Prentice Hall; 1976.

[16] Timoshenko S, Young DH. Vibration problems in engineering. Princeton, NJ: D. Van Nostrand Company, Inc.; 1955.

[17] Lee Y, Feng ZC. Dynamic responses to sinusoidal excitations of beams with frictional joints. Commun Nonlinear Sci Numer Simul 2004;9:571–81.