On the general theory of chaotic dynamics of flexible curvilinear Euler–Bernoulli beams

We present chaotic dynamics of flexible curvilinear shallow Euler–Bernoulli beams. The continuous problem is reduced to the Cauchy problem by the finite-difference method of the second-order accuracy and finite element method (FEM). The Cauchy problem is solved through the fourth- and sixth-order Runge–Kutta methods with respect to time. This preserves reliability of the obtained results. Nonlinear dynamics is investigated with the help of a qualitative theory of differential equations. Frequency power spectra using fast Fourier transform, phase and modal portraits, autocorrelation functions, spatiotemporal dynamics of the beam, 2D and 3D Morlet wavelets, and Poincaré sections are constructed. Four first Lyapunov exponents are estimated using the Wolf algorithm. Transitions from regular to chaotic dynamics are detected, illustrated and discussed. Depending on signs of four Lyapunov exponents the chaotic, hyper chaotic, hyper-hyper chaotic, and deep chaotic dynamics is reported. Curvilinear beams are treated as systems with an infinite number of degrees of freedom. Charts of vibration character, elastic–plastic deformations, and stability loss zone versus control parameters of the studied beams are reported.


Introduction
Nonlinear vibrations of Euler-Bernoulli beams have been studied for a long time. An Euler-Bernoulli beam represents a continuous structural member and its vibrations are governed by nonlinear partial differential equations, and hence, there is no hope to find their exact analytical solutions. Therefore, the solutions have been detected in an approximate way using either analytical (mainly asymptotic) or numerical methods. In this paper, we solve initial boundary-value problems using numerical techniques (a discussion on the advantages and disadvantages of numerical versus analytical approaches can be found elsewhere [1][2][3]).
It is clear that various types of beams, i.e. linear and nonlinear, uniform and non-uniform, prismatic and non-prismatic, etc. have been modeled and analyzed intensively in many mechanical, civil engineering and mechatronic structures/mechanisms/devices such as buildings, bridges, ships, turbines, helicopter blades, aircraft measurement and control devices, sensors and transducers just to mention a few of them.
There is a large number of papers devoted to the analysis of beam vibrations, and also a rather wide spectrum of different directions of research dedicated to these problems. The reason is motivated mainly by relatively simple PDEs that govern the dynamics of beams, and hence, various numerical analytical methods can be applied in order to achieve reliable and validated results.
On the other hand, various types of beam shapes and design can be taken into account along with various types of nonlinearity, including geometric, physical and design.
Below, we briefly describe the state of the art of recent directions in modeling and analysis of beam vibrations putting emphasis on a novel research direction proposed in our paper.
The Euler-Bernoulli beam with simply supported fixed ends modeling power transmission belts was studied both theoretically and experimentally by Moon and Wickert [4]. Pitchfork bifurcations, the stability of equilibria, as well as post-bifurcation velocity range were detected, illustrated and discussed. This study was extended by Pellicano and Vestroni [5], where the rub and super-critical speed intervals and their influence on the stability, bifurcation and global nonlinear dynamics of the Euler-Bernoulli axially moving beam were extensively analyzed. Exact solutions in closed form of forced vibrations of prismatic damped Euler-Bernoulli beams as well as Green functions for the beams with different homogeneous and elastic boundary conditions were reported by Abu-Hilal [6]. Naguleswaran [7] investigated the transverse vibration of a homogeneous Euler-Bernoulli beam of constant depth and linearly varying breadth taking into account damped, pinned, sliding and free boundary conditions, and derived associated frequency equations. Park and Gao [8] proposed a novel model for the bending of an Euler-Bernoulli beam using a modified couple stress theory.
A new model for prediction of bending rigidity of the cantilever beam versus the classical beam model was reported. Bifurcations and orbital stability loss exhibited by N nonlinearly coupled Euler-Bernoulli microelectromechanical beams parametrically actuated were studied by Gutschmidt and Gottlieb [9]. They compared the analytically obtained results with numerical ones and complemented them with a three-beam array analysis. A new semi-analytical method to study nonlinear vibrations of Euler-Bernoulli beams with general boundary condition was proposed and tested in reference [10]. Recently, both variational iteration and parameterized methods have been applied to study nonlinear vibration of Euler-Bernoulli beams subjected to axial loads [11]. This type of study has been extended to analyze nonlinear responses of a damped-damped buckled beam [12].
The Euler-Bernoulli curved beams are less investigated. Lim et al. [13] derived exact relationships between deflections and stress resultants of Timoshenko and Euler-Bernoulli curved beams. Hiller [14] analyzed three elastic degrees of freedom for each link as an Euler-Bernoulli beam assuming inner constraints for the other three elastic co-ordinates. Shi et al. [15] derived a complete second-order deformation field together with the governing equations within a graphtheoretic formulation of the beam model. An enhanced nonlinear 3D Euler-Bernoulli beam model governed by ten coupled PDEs was proposed by Zohoor and Khorsandijou [16], and then, it was extended with newly added elastic terms exhibiting tension-compression, torsion and two spatial bending dynamical effects [17].
Our paper addresses a novel approach from the point of view of reliable and validated results of strongly nonlinear governing PDEs and also in reference to modern approaches to the detection and monitoring of regular and chaotic dynamics of this type of structural members. It extends our earlier studies reported in references [18][19][20][21][22][23][24].
The paper is organized in the following manner. In Sect. 2, equations governing nonlinear dynamics of curvilinear Euler-Bernoulli beams are given in the non-dimensional form and the initial boundary-value problem is defined. Quantitative and qualitative analyses of the results obtained via finite-difference method (FDM) and finite element method (FEM) are carried out in Sects. 3 and 4, respectively. Section 5 is devoted to the computation of four Lyapunov exponents, whereas the last Sect. 6 summaries the obtained results. In this work, we consider flexible one layer thin curvilinear beams of length a, height h and geometric curva- where R x is the beam curvature radius. The beam is loaded continuously along its external surface q(x, t), acting in the normal direction to its middle surface ( Fig. 1).
A mathematical model of the investigated beam is based on the following hypotheses: (i) An arbitrary transversal beam cross section normal to its middle surface is straight and normal before and after a deformation, and hence the beam cross section height remains unchanged; (ii) Although the inertia of rotation of beam elements is neglected, inertial forces, being responsible for displacements along both normal and tangential directions, are taken into account; (iii) External forces do not change their directions during beam deformations; (iv) Geometric nonlinearity is taken in Kármán's form; (v) Thin curvilinear beam exhibits shallow properties.
The hypotheses listed so far are based on the Euler-Bernoulli ideas, and though they yield a mathematical model of the first-order approximation, it is sufficiently rigorous for the analysis of thin curvilinear beams.
A mathematical model of the beam is governed by a system of nonlinear PDEs. They describe the way of beam element movement taking into account energy dissipation, and they are given in the following nondimensional form [25] where L 1 (u, w), L 2 (w, w), L 3 (w, w) are the nonlinear operators; w(x, t) is the beam normal deflection; u(x, t) is the beam element longitudinal displacement; ε 1 -coefficient of dissipation of a surrounding medium; E-Young modulus; h-height of the transversal beam cross section; γ -unit weight density; g-Earth acceleration; k x -beam curvature; ttime; q = q 0 sin(ω p t)-external load; q 0 -amplitude of the external load; ω p -external load frequency. Non-dimensional parameters are introduced in the following way: and bars over non-dimensional quantities are omitted in Eq. (2.1). The following boundary condition is supplemented to Eq. (2.1): which means that one bar end is fixed (x = 0), whereas the second one is pinned (x = a), and the following initial conditions are taken: In what follows we report numerical data for boundary conditions (2.3) and initial conditions (2.4) taking into account periodic excitation q = q 0 sin(ω p t). Equations (2.1)-(2.4) cannot be solved analytically, and hence, they will be solved numerically applying the FDM of the second-order accuracy with respect to spatial coordinate x [26]. In order to reduce PDEs to ODEs with respect to time co-ordinate, we use the finite-difference approximations applying Taylor series in the neighborhood of point x i . We consider a mesh The following difference operators with the approximation of O(c 2 ), where c is the step of spatial coordinate partition, are introduced: Hence, PDEs (2.1) are reduced to the second-order ODEs regarding time co-ordinate in the form After the application of the second-order finitedifference approximation, the obtained ODEs (2.5) with the corresponding boundary and initial conditions are solved using either the fourth-or sixth-order Runge-Kutta method [27]. Reliability and validity of the obtained results have been followed via two qualitatively different approaches, i.e. FDM and FEM. A comparison of signals (time histories), Fourier power spectra, and Morlet wavelets proved convergence of the results obtained by both FDM and FEM. It should be noted that we use further a standard fourth-order Runge-Kutta method to solve Eq. (2.5), since the sixth-order Runge-Kutta scheme yielded the same results while requiring twice as long computational time.
In order to investigate the dynamics of a thin curvilinear flexible beam, the algorithms and associated programs were developed to construct and monitor signals, frequency power spectra, phase and modal portraits, Poincaré sections, autocorrelation function and the Lyapunov exponents. Using the developed programs, it is possible to choose q 0 , ω p to define and monitor various types of vibrations (periodic, period doubling bifurcations, quasi-periodic and chaotic vibrations, etc.). Graphic interpretation of the results is given in charts exhibiting beam vibration types versus control parameters q 0 , ω p .

Numerical analysis via charts
We consider a curvilinear beam with the following parameters n = 120, λ = 100, k x = 48 and ε 1 = 1. The constructed charts of vibration regimes not only exhibit a picture of a studied nonlinear process, but allow us to recognize validity and importance of the chosen intervals of particular vibration zones. In this case, they present the dependence of a vibration regime on control parameters, here frequency ω p and excitation amplitude q 0 . An optimal choice of the number of nodes n applied in the FDM mesh, the number of partition of the investigated interval of excitation frequency ω p , and the excitation amplitude is equal to 90. Figure 2 shows charts and color notation regarding the type of curvilinear beam vibrations and zones of qualitatively different vibrations with respect to control parameters {ω p , q 0 }.
Besides the type of beam vibrations, we are also interested in a series of quantitative characteristics of the beam vibration process. This important information is provided by bounded beam deflections introduced via the hypotheses while deriving a mathematical model of the beam. This is why construction of the chart requires monitoring of maximum (absolute value) deflection. Figure 3 shows two charts of threshold regimes, i.e. those allowed and not allowed owing to the applied hypotheses. The first chart includes absolute deflections, whereas the second one presents a derivative of the first one.
It should be emphasized that for large zones the obtained deflections are out of the assumed hypotheses (the chart in Fig. 3a has a rather theoretical meaning and explains character of the chart reported in Fig. 3b). In Fig. 3, black (white) color corresponds to a minimum Another valid characteristic is a threshold extension of the curvilinear beam. In order to construct charts of threshold extensions, each numerical experiment is supplemented with the maximum extension of the curvilinear beam. The latter characteristic allows us to detect zones of plastic-elastic deformations. Introduction of relative beam extension ε = l 1 −l l , implies Hook's law in relative units σ = Eε, where E is the elasticity modulus. Knowing the limit of proportionality and offset yield strength for chosen materials one may construct charts with zones of elasticplastic deformations. Figures 4 and 5 present charts of threshold deformations and a chart with elastic-plastic zones for thermally strengthened aluminum 1915T (this material is chosen owing to its high upper yield point).
Averaged deflection per period yields a displacement of the middle beam curve. Figure 6 shows the chart of averaged deflections, where red (blue) color corresponds to positive (negative) maximal beam deflection, and the black one reflects the lack of dis- placement. Owing to the reported chart, a rapid displacement of the averaged beam surface on the defined boundary occurs, which implies dynamical stability loss of the beam. Another important result relating to nonlinear analysis is yielded by a chart of signal divergence which is constructed in the following way. Let a i and b i be the deflections regarding "a" and "b" series, respectively. Hence, one may define a distance between signals using the formula where M stands for the number of points in a signal. As a result, we get a scalar quantity characterizing divergence (difference) between the series, being evaluated through a metric applied to the set of signals with one frequency of excitation. In order to avoid an abstract quantity measuring the distance between series, we transit to formula S = ( M i=1 |a i − b i |)/M. Therefore, the averaged distance (quantified in deflections) between points of the series "a" and "b" is obtained. If the distance tends to zero, then signal "a" tends to signal "b", and if the distance is equal to zero then we deal with one signal only. If a curvilinear beam (system of equations) exhibits a linear behavior, then a uniform increase of the load yields uniform divergence of signals. On the contrary, if a linear behavior is not observed, then a space of nonlinear behavior is fixed. This kind of dependence is constructed for a load acting on the beam. The latter is responsible for the amount of energy entering the system. Construction of this dependence is rather difficult with respect to frequency of the exciting load owing to two main reasons: (1) the frequency is responsible for both the observed vibration regime and the way of energy transmission into the system; (2) the frequency appears in a formula governing the excitation load which is usually described by a trigonometric function, i.e. an essentially nonlinear relation. Therefore, construction of the vibration chart requires transition along a vertical line (up and down), and the monitored signal is compared with the previous one in order to fix the estimated divergence. Figure 7 shows charts of signal divergence in absolute units and with maximum cuts of the amount of 2, 0.5, and 0.25 deflection units (white color-signal divergence exceeded the assumed maximum, green color-result of modeling not defined). Figure 8 illustrates the chart of vibration regimes and the chart of signal divergence with a maximum cut of 0.25 units.
On the basis of data reported in Fig. 8, the following conclusions can be formulated: 1. Charts of vibration regimes and of signal divergence are self-interacting and self-completed. 2. Series of hypotheses are applied while establishing the beam mathematical model. Hence, putting limitations to threshold deflections of the beam, one may define the space of validity and reliability of the applied mathematical model. Analysis of the results beyond the space defined in this way yields also a theoretical contribution to the studied problem. 3. The chart of threshold extension defines zones of elastic-plastic deformations for different materials. In Fig. 8, we deal with the thermally reinforced aluminum 1915T. It is evident that the zone of elastic deformations almost coincides with the space of validity of the applied mathematical model. For comparison, Fig. 9 shows charts of vibration regimes for the following fixed parameters n = 120, λ = 100, ε 1 = 1, and for k x = 0, 12, 48 (k x = 0 corresponds to a classical beam). One may conclude that the beam curvature increase implies an increase of zones with periodic vibrations. The increase of beam curvature causes an increase of beam stiffness which yields an increase of the zones of periodic vibrations. Figure 10 shows charts of threshold deflection values (cut of 5 units).
It can be concluded that an increase of the beam curvature implies a decrease of applicability of the used mathematical model. Besides, a boundary (threshold) of applicability is clearly visible through a sudden change of the vibrations regime, i.e. dynamical stability loss occurs. Figure 11 includes charts of signal diver-  One may conclude that an increase of the beam curvature implies the increase of zones with dynamical stability loss.

Further numerical investigations
A beam with the curvature of k x = 48 and parameters n = 120, λ = 100, ε 1 = 1 is further studied. Analysis is carried out for the fixed excitation frequency value (ω p = 5, 7615) and various external load amplitudes. In each component of Fig. 12 It is now well recognized that a wavelet transform allows us to investigate time-frequency dependencies of various dynamical systems. In this paper, we apply the Morlet wavelet transform. More descriptions of the application of different wavelets to study spatiotemporal dynamics of structural members can be found in references [19,28,29], where advantages of the use of Morlet wavelet are illustrated and discussed in particular.
Analysis of the results given in Fig. 12 leads to the following conclusions. For the excitation amplitude q 0 = 57,500 (see Fig. 12a), the beam exhibits periodic vibrations. The scale of the attractor is so small that it does not influence the Fourier power spectrum. A relatively small increase of the excitation amplitude up to q 0 = 62,500 (see Fig. 12b) generates one independent and a set of dependent frequencies; phase portrait exhibits an attractor, the modal portrait is collapsed, and the Poincaré section exhibits the attractor. For the amplitude interval q 0 = 64,500 (see Fig. 12c)-83,000 the situation does not change qualitatively, though slight changes take place in the Poincaré section as well as in the energy distribution along frequencies in the power spectrum. For q 0 = 85,000 (see Fig. 12d), a collapse of both previous Poincaré section and phase portrait is observed, and a step-wise transition into a periodic vibration regime is manifested for q 0 = 86,000. However, this state is not robust, and already in the interval of q 0 = 86,500 (see Fig. 12e)-87,000 (see Fig. 12f) a transition into chaotic dynamics is observed (2D Morlet wavelet exhibits an intermittency of frequencies, whereas the 3D Morlet wavelet shows a transition of the energy into the zone of low frequencies; Poincaré section transits into a new attractor, phase and modal portraits exhibit this attractor, and the autocorrelation function rapidly decreases). It should be emphasized that changes discussed so far refer not only to the vibration character but also to such important characteristics as the maximum deflection. Observe that for the excitation amplitude q 0 = 86,500 (see Fig. 12e) it is equal to 1.75, for q 0 = 87,000 (see Fig. 12f) it is 4.25, whereas for q 0 = 88,500 (see Fig. 12g) it is equal to 10.5. The maximum deflection moves to the beam center along the co-ordinate x, and the system loses its sensitivity to non-symmetric boundary conditions. In practice, our investigated thin curvilinear beam exhibits buckling and loses a possibility to carry the external load, which is well illustrated through the reported beam deflection.
A further analysis is mainly theoretical/speculative one, because the deflection of the curvilinear beam exceeded the value assumed in the applied hypothe-  Analysis of the Fourier frequency power spectrum implies the following scenarios of transition into chaos. For the excitation amplitude lesser than q 0 = 57,500 periodic orbits with frequency ω p = 5.7615 are observed. An increase of the excitation amplitude causes the Sacker-Neimark bifurcation, and frequency ω 1 = 2.614 non-commensurable with the fundamental one appears. In practice, the simultaneous birth of linearly dependent frequency ω 2 = ω p − ω 1 is observed. Next, the Ruelle-Takens-Newhouse scenario [30] is exhibited, which is clearly demonstrated for q 0 = 62,500 (see Fig. 12b). This scenario takes place up to q 0 = 86,000, and then the beam dynamics becomes periodic with only one frequency (har-monic). Increasing the load from q 0 = 86,500 (see Fig. 12e) up to q 0 = 87,000 (see Fig. 12f) an intermittency phenomenon appears (this is well reported by the 2D Morlet wavelet), which yields a noisy frequency spectrum. In spite of that three frequencies ω 1 = 2.405, ω p , ω 2 = ω p − ω 1 are distinguishable, and again the Ruelle-Takens-Newhouse scenario occurs. This classical version of the Ruelle-Takens-Newhouse scenario takes place up to q 0 = 87,500. A further increase of the excitation amplitude implies an interesting beam dynamical behavior, i.e. it exhibits simultaneously two scenarios. Namely, for q 0 = 88,500 (see Fig. 12g) we may monitor fundamental frequency ω p , two frequencies ω 1 and ω 2 = ω p − ω 1 from the previous Ruelle-Takens-Newhouse scenario, and frequencies ω 3 = ω p /2 = 2.884, ω 4 = ω p /4 = 1.424, ω 5 = 3/4 ω p = 4.308. Since the values of frequencies ω 3 , ω 4 , ω 5 are obtained through Hopf bifurcations, Feigenbaum's scenario (see, for example, Chapter 9 of monograph [31]) governs the beam dynamics while approaching chaos. Then, with a further increase of excitation amplitude, the system prolongs this two-state dynamics. However, for q 0 = 95,500 (see Fig. 12j), the beam dynamics is qualitatively changed, since the following frequencies, as well as frequencies ω 2 = 1.387, ω 3 = 1.497, ω 4 = 4.271 and ω 5 = 4.369 appear. It is remarkable that frequencies ω 2 and ω 3 are concentrated around ω p /4 = 1.436, whereas frequencies ω 4 and ω 5 around 3/4ω p = 4.307. In other words, the Hopf bifurcation in the implicit and explicit forms, as well as the damped Ruelle-Takens-Newhouse scenario are observed. A more detailed analysis of the beam vibrations within the interval from q 0 = 95,000 (see Fig. 12h) to q 0 = 96,000 showed that the beam moved between two illustrated dynamical states a few times in a step-wise manner. Two-state system dynamics is robust and also appears for higher values of q 0 , until it reaches a deep chaotic state with a broad band of Fourier spectrum.

Lyapunov exponents
As it is well known, the Lyapunov exponents are used to estimate quantitatively the order of dynamical chaos of a studied system. The number of Lyapunov exponents depends on the number of differential equations governing the system dynamics. Each equation is associated with one Lyapunov exponent. If a Lyapunov exponent is negative, then we deal with a periodic process, and when we have a positive exponent, then a chaotic component within the system trajectory occurs. Although there are different methods to compute the Lyapunov exponents, we apply the Wolf method here [32]. This method allows us to compute the largest Lyapunov exponents by choosing only one coordinate (time history), particularly suitable when the system of governing equations is not known and we cannot measure all of their coordinates. Let us have a time series x (t) , t = 1, . . . , N yielded by the measurement of one coordinate of the chaotic process, and the successive measurement carried out for equal time distances. Hence, the method of mutual information is used to define time delay τ , whereas the method of false neighborhoods yields the estimation of an embedding space dimension m. The so far described reconstruction process yields a set of points of space R m : (5.1) where i = ((m − 1) τ + 1) , . . . , N . Let us take a point from series (5.1) and denote it by x 0 . Monitoring of series (5.1) allows us to find such point x 0 that the following inequality holds x 0 − x 0 = ε 0 < ε, where ε is the fixed quantity, much lesser than the dimension of a reconstructed attractor. It is necessary that points x 0 and x 0 be separated in time. Then, evolution of these two points is followed until the distance between them achieves an a priori given quantity ε max . We denote the obtained points by x 1 and x 1 , a distance between them by ε 0 , and the interval of time evolution by T 1 . After that we again consider series (5.1), and we find point x 1 close to x 1 , requiring the following non-equality to be satisfied and x 1 − x 1 should possibly have the same direction. Further, the procedure is extended to points x 1 and x 1 . Repeating the so far described procedure a few times M, the largest Lyapunov exponent is estimated through the following formula A more detailed description of Wolf algorithm can be found in reference [32]. In this work, the application of Wolf's method yielded Lyapunov exponents for signals monitored in the beam middle point and for different values of the excitation load amplitude q 0 (we have investigated interval of q 0 ∈ [500; 199,500] with the step of 500 units). The obtained results are reported in Fig. 13. In addition, a computation of the four Lyapunov exponents was carried out for the signals given already in the above. Figure 13 illustrates a transition of the beam through a periodic zone (all exponents are either less or equal to zero), hybrid zone of system states (third and fourth exponents change their values keeping mainly positive values), and then it transits into a deep chaotic state with all four exponents positive. These results coincide in full with the analysis carried out in the previous part of the paper. It should be mentioned that in spite of the described way of analysis we have also constructed the so-called charts of the Lyapunov exponents. The charts of Lyapunov exponents, similarly to the charts illustrated in Sect. 3, are constructed versus control parameters ω p , q 0 with the same intervals as in the mentioned section. Figure 14 shows three charts of the Lyapunov exponents.
Comparison of the charts of Figs. 2 and 14c shows qualitatively similar results perhaps in spite of the estimation of quasi-periodic regimes. However, it should be noted that quasi-periodicity detection using the Lyapunov exponents is less sure in comparison to other characteristics applied while constructing the chart of Fig. 2, hence the latter one is more reliable.
Comparison of the results shown in Fig. 14 with the results included in Figs. 8 and 9 implies high coincidence obtained via the qualitatively different ways of the analysis carried out. Figures 12 and 13 confirm and complete the predictions from Fig. 2. The case of presented windows of the q 0 choice associated with quasi-periodic regimes shown in Fig. 12 corresponds to the case when all of the four Lyapunov exponents are equal to zero. Furthermore, each chart form shows a different aspect of the dynamical process, and matching the results obtained in different ways give not only a reliable and validated picture of the system behavior but also helps to understand it fully.

Concluding remarks
The flexible curvilinear beams are widely applied in robotic manipulator arms, helicopter rotor blades, aircraft propellers, flexible satellites, multi-packet blade systems and wind turbines. The general theory of flexible beam nonlinear dynamics proposed in this paper gives a recipe for direct engineering applications. Namely, the constructed charts of the beam vibration regimes allow us to avoid dangerous ones, and hence give hints to control regular and chaotic dynamics of the studied beams.
In spite of the mentioned impact of the results on the applications of flexible Euler-Bernoulli beam, the paper offers a novel insight into reduction of the derived partial differential equations to a set of nonlinear ordinary differential equations putting emphasis on the reliability and validation of the obtained results. We have applied not widely used approaches to study bifurcation and chaotic beam vibrations by 2D and 3D Morlet wavelets and four Lyapunov exponents computed via the Wolf algorithm along with the classical numerical tools like frequency power spectra, phase and modal portraits, autocorrelation functions and Poincaré maps.
In addition, the paper includes illustration and discussion of novel physical nonlinear phenomena yielded by the studied flexible beams, which can be exhibited also by other structural members. The research fits also plastic-elastic beam deformations. Since the paper is written in a concise style, we omit here repeating of the whole detected and illustrated changes of the vibrational beam regimes implied by the changes of the chosen control parameters, i.e. amplitude and frequency of the continuously loaded harmonic excitation.
Let us emphasize only that our multi degree-offreedom nonlinear mechanical systems exhibit either standard scenarios of a route from periodicity to chaos (Feigenbaum's and the Ruelle-Takens-Newhouse scenarios), or a simultaneous occurrence of both mentioned classical scenarios while entering the full chaotic regime. Finally, it is worth mentioning that detection of chaotic beam regimes associated with the so-called hyper-hyper chaos and deep chaos (three and four Lyapunov exponents are positive, respectively) belongs certainly to new results and its validation via a laboratory experiment is highly demanded.