Chaotic vibrations of flexible shallow axially symmetric shells

In this work, chaotic dynamics of flexible spherical axially symmetric shallow shells subjected to sinusoidal transverse load is studied with emphasis put on the vibration modes. Chaos reliability is verified and validated by solving the implemented mathematical model by partial nonlinear equations governing the dynamics of flexible spherical shells and by estimating the signs of the largest Lyapunov exponents with the help of qualitatively different approaches. It is shown how the scenario of transition of the investigated shells from regular to chaotic vibrations depends on the boundary condition. The following cases are considered: (1) movable and fixed simple supports along the shell contours, taking into account shell stiffness (Feigenbaum scenario) and shell damping (Ruelle–Takens–Newhouse scenario), and (2) movable clamping (regular shell vibrations). The presence of dents, the location and character of which essentially depend on the shell geometric parameters, boundary conditions, and the external load parameters, is detected in some regions of the shell surface and discussed.


Introduction
Axially symmetric spherical shells, which are examples of thin-walled constructions, are widely applied in aviation and rocket industries, shipbuilding, automotive industries, energy-harvesting and chemical industries, fabrication of sensor devices, and bioengineering. In numerous constructions, the mentioned structural members are subjected to different external loads and boundary conditions, and thus, it is important to study their nonlinear vibrations.
The dynamics of flexible shells has a long history in mechanics, and only a selected part of it is considered in the present paper, putting emphasis on the novelty of the present work with respect to the findings in the available literature.
Baker [1] presented a theory predicting the existence of two infinite sets of the proposed normal shell model having bounded/unbounded frequencies for studying the axisymmetric modes of vibrations of thin spherical shells.
The transverse vibration of shallow spherical shells was investigated by Kalninis and Naghdi [2] and Kalninis [3,4], including the shear theory and free vibration problems. Jain [5] analyzed the axisymmetric vibration of a shallow spherical shell subjected to timevarying loads, including the effects of transverse shear and rotary inertia. Special cases of the concentrated ring and central loads as well as the symmetrically distributed loads were also addressed.
The investigations devoted to studying the character of dents exhibited by shell surfaces were carried out by Volmir [6], Pogorelov [7], and Grigoluk and Kabanov [8], whereas the dents localization was analyzed by Mikhasev and Tovstin [9] with the help of the asymptotic methods.
Okazaki et al. [10] carried out a theoretical analysis of damping of the flexural vibration of a threelayer shallow spherical shell with a constrained viscoelastic layer. The influence of changes in the opening angle, shear parameter, layer thickness, and introduced boundary conditions on shell vibration was studied in terms of Bessel functions.
Evkin and Kalamkarov [11,12] employed a perturbation analysis to large deflection of equilibrium states of composite shells of revolution. Evkin [13] proposed an asymptotic solution to the problem of large deflections of deep orthotropic spherical shells subjected to radial concentrated load.
Creep buckling of shallow arches was studied by Wang et al. [14].
Hamed et al. [15,16] investigated creep buckling of imperfect thin-walled shallow concrete domes and nonlinear long-term behavior of spherical shallow shells of revolution.
Li et al. [17] analyzed the effect of spherical shell geometry and material properties on the critical normal load, interference, and the contact area where plastic yielding occurred. The authors proposed an empirically found relation for the load-interference spherical shell behavior prior to its plastic yielding. Gurijala and Perati [18] studied axially symmetric vibrations of composite poroelastic spherical shell composed of two spherical shells. Biot's theory was employed to derive the frequency equations for pervious/impervious shell surfaces. The dependence of nondimensional phase velocity on the ratio of outer and inner radii was estimated for two types of sandstone spherical shells.
Shi et al. [19] detected and studied the effects of curvature on the vibration characteristics of doubly curved shallow shells with general elastic edge constraints. A unified novel solution was proposed and validated by numerous tabular and graphical results illustrating the influence of curvature on the shell material frequencies for different boundary conditions. However, no studies aimed at analyzing the modes of vibration of flexible circle axially symmetric shells subjected to uniformly disturbed harmonic loads have been found in the available, published literature. Furthermore, to the knowledge of the authors, the occurrence of localized dents exhibited by the spherical axially symmetric shells and their matching with the chaos theory is reported for the first time in the present manuscript.
Furthermore, novel scenarios of transition from regular to chaotic vibrations of the shell under the sinusoidal transverse load are presented for four boundary conditions (simple movable/nonmovable supports and movable/nonmovable clamping).
Our investigations are aimed at detecting "true" chaotic vibration. The vibration forms as well as the occurrence of located dents exhibited by flexible axially symmetric shallow shells depending on the employed boundary condition are studied.
The paper is organized in the following way: Sect. 2 deals with the problem statement, the mathematical model, and reduction from PDEs to ODEs by means of the FDM (finite difference method). Section 3 is devoted to studying the reliability and validity of the shell chaotic vibrations, whereas modes of nonlinear shell vibrations are investigated in Sects. 4 and 5.

Problem statement and the method of solution
Let a spherical shallow shell be considered in the polar system of coordinates introduced in the following way: The studied axially symmetric shallow shell with the uniformly distributed harmonic load (a) and the scheme of a thin shell (b) (see Fig. 1a). Owing to the earlier work of Novozhilov [20], a shell is called thin if h/R 1 (see Fig. 1b), and the inequality h/R ≤ 1/40 should be employed in such a case, as suggested by Novozhilov. If this inequality is not satisfied, then the shell is called thick. However, the majority of practically applied shells satisfy the inequalities 1/2000 ≤ h/R ≤ 1/50, i.e., they belong to the class of thin shells. On the other hand, the threshold used to identify shallow shells is defined by the Reissner [21] ratio f /2c << 1/8, where f is the distance between the shell vertex and a plane f and the distance 2c is measured in plane (see Fig. 1b). Vlasov [22] suggested f /2c ≤ 1/5. Employing the Reissner criterion for shallowness (α < α * = 60 • ), one gets 2c ≤ R.
The system of nonlinear PDEs governing dynamics of the axially symmetric shells has the following form [23]: The below nondimensional quantities (with bars) are introduced:  [24,25] and Appendix A). In the present investigations, we are aimed at studying the boundary conditions and two excitation parameters q 0 and ω p of the sinusoidally distributed q = q 0 sin ω p t.
The following boundary conditions are considered: (1) simple contour movable in the meridional direction In order to reduce problem (2.1)-(2.7) governing dynamics of the considered continuous system to a system with lumped parameters, the finite difference method (FDM) with approximation O( 2 ) is used. PDEs as well as boundary and initial conditions (2.1)-(2.7) are recast to the following finite difference formulas with respect to the spatial coordinate r and time: where = b/n and n denotes the number of nodes of the shell radius.
The finite difference method (FDM) is employed only to spatial coordinates. It reduces the problem governed by PDEs (2.1) of an infinite dimension to the system of ODEs (2.8) with respect to time (the Cauchy problem). The Cauchy problem is solved using a few methods of the Runge-Kutta types, i.e., the fourth-order Runge-Kutta method (RK4), the secondorder Runge-Kutta method (RK2), the fourth-order Runge-Kutta-Fehlberg method (RKF45), the fourthorder Cash-Karp method (RKCK), the eighth-order Runge-Kutta Prince-Dormand method (RK8pd), and the second-order and the fourth-order implicit Runge-Kutta method (rk2imp and rk4imp). After conducting numerous case studies, we have finally chosen the fourth-order Runge-Kutta method. The reason is that the accuracy of the this method is the same as of the Runge-Kutta methods of higher orders, but the computational time is shorter. Location of nodes along the shell radius r is shown in Fig. 2.
The counterpart difference forms of the boundary conditions are as follows: 1. Simple contour movable in the meridional direction 2. Rigidly clamped contour (2.10) 3. Sliding clamping of the shell contour (2.11)

Simple nonmovable contour
= w n−1 , w n = 0 for r n = b; (2.12) and the below initial conditions are taken (2.13) If small terms are neglected and differential operators are substituted by the central finite differences for r = , the following conditions are obtained at the top of the shell: 14) The FDM (the finite difference method) is an alternative to the Faedo-Galerkin, the Ritz, and the finite element methods. Let us cite a few recently published papers, in which the FDM has been employed to study the dynamics of structural members. Medina et al. [26] studied axially symmetric nonlinear behavior of electrostatically actuated and initially curved bistable microplates. Chen et al. [27] considered post-buckling of size-dependent microplate, focusing on damage effects due to the transverse load. Mehditabar et al. [28] carried out the thermoelastic study of a functionally graded piezoelectric rotating hollow cylindrical shell subjected to dynamic loads.
Aranda-Iglesias et al. [29] also investigated spherical shells. The solution was found using both the finite element method and the finite difference method. The same results were obtained for both methods, which validates the reliability of the FDM.
The transverse load can be changed arbitrarily with respect to the spatial coordinate and time. In this work, the harmonic transverse load of the form q = q 0 sin(ω p t), where q 0 stands for an amplitude and ω p is a frequency of the excitation, is used.

Quantification of true chaotic vibrations
The problem related to reliability of chaotic vibration is extremely important, since in many cases errors intro-duced by numerical computations can be interpreted as chaotic behavior, which has been illustrated and discussed in reference [30]. It is known that one of the fundamental properties of chaotic dynamics is associated with its essential dependence on the initial conditions. On the other hand, more rigorous definition of chaos, given by Duavaney in 1989 [31], aside from the already mentioned property, emphasizes the importance of two other characteristics of chaos, i.e., intermittency/transitivity and regularity/periodicity conditions manifested by the density of periodic points. In 1992, Banks et al. [32] proved that the transitivity and periodicity conditions imply sensitivity to the initial conditions. On the other hand, Knudsen [33] defined chaos with the help of a function given on the bounded metric space, which is understood as a chaotic one if it exhibits a dense orbit and presents an essential dependence on the initial conditions.
Owing to Gulick results [34], chaos exists if sensitivity to the initial conditions is exhibited or a positive LE (Lyapunov exponent) is identified in each point of the space on which the chaotic function is defined, and the latter does not tend to periodicity. In this work, we employ the Gulick definition.
To solve the problem of reliability of chaotic vibrations while numerically solving nonlinear PDEs governing dynamics of spherical axially symmetric shells subjected to harmonic transverse load, the following steps are carried out: 1. Since PDEs are reduced to a Cauchy problem by means of the second-order FDM, the solution strongly depends on the number of nodes involved in the shell radius r ∈ [0; b] partition. In the present paper, the solution is accepted if for n and n + 1 nodes the obtained solutions coincide. In particular, the convergence of the employed FDM is analyzed for different shallowness parameters and boundary conditions. 2. Since the Cauchy problem is also solved numerically, the solution essentially depends on both the used method and the employed time step. Therefore, the following methods are used: fourthand second-order Runge-Kutta (RK4 and RK2) methods, the fourth-order Runge-Kutta-Fehlberg (RKF45) method, the fourth-order Cash-Karp (RKCK) method, the combination of the eighthorder Runge-Kutta and Dormand method (RK8pd), as well as the second-(rk2imp) and the fourth-order (rk4imp) implicit Runge-Kutta methods. 3. The numerically obtained data are used to construct the time series (signals), phase portraits, Fourier power spectra (FFT), Morlet wavelets, and shell deflection surfaces. The following wavelets are used: Meer wavelets, Shannon-Kotelnikov and Meier wavelets, Daubechies wavelets from db2 to db4, Coiflets, Morlet wavelets, complex Morlet wavelets as well as the wavelets based on the derivative of the Gauss function of order higher than 8. It should be emphasized that it was not feasible to investigate the considered problem with Meer and Shannon-Kotelnikov wavelets, i.e., the first one did not identify the frequency zones properly, whereas the second one did not yield the expected localization in time. The carried out analysis of the wavelet spectra obtained by the Daubechies wavelets, Coiflets, and simlets has shown that an increase in the order of the employed filter improves proper frequency identification. Furthermore, in spite of the differences between the employed wavelets and their filters, the wavelet spectra yielded by the Daubechies wavelets, simlets, and Coiflets of the same order give almost the same results. However, the frequency localization yielded by them is not accurate enough, which makes them unsuitable for analysis of vibration of the investigated mechanical system. Considering the results obtained by means of derivatives of the Gauss functions, the accuracy of the frequency localizations increased with the increase in the order of the derivative. It should be noted that the frequency spectra obtained by the Meer wavelet, which is a smooth variant of the Shannon-Kotelnikov wavelet, have been localized better than in the case of Morlet wavelets only for a low-frequency band. In what follows we give results for Morlet wavelets that were specially chosen and adopted to study the stated problem. In these wavelets, there is no scale function ϕ, and ψ is given explicitly, i.e., there is no compact carrier. Therefore, to study complex vibration of continuous mechanical system, like a plate/shell, one can employ either complex or real Morlet wavelets or wavelets based on either real or complex derivatives of the Gauss functions of the order higher than 16. 4. Since the definition of chaos given by Gulick [34] is used in this work, one needs to compute Lya-punov exponents (LEs) and define their signs. In this paper, a spectrum of the LEs for each partition of the interval [0; b] was obtained by three methods: Wolf's [35], Rosenstein's [36], and Kantz's [37]. Furthermore, convergence of the above-mentioned methods was also investigated.
The final solution was accepted and treated as reliable/true if coincidence up to four digits after the decimal point was achieved. According to the already given steps of investigation, convergence of the solution depending on the number n of partitions of the shell radius was considered at first. The study was carried out for all types of boundary conditions: movable simple support (2.9), rigid clamping (2.10), movable clamping (2.11), and nonmovable simple support (2.12) for b ∈ [3; 8]. Then, results for n = 10, 20, 40, 80, 100, 120, 140 were compared. After reducing problem (2.1)-(2.5) to its counterpart normal form, the first equation of the Cauchy problem was solved by the Runge-Kutta methods with respect to the deflection function w. At the same time, algebraic equations were solved (on each time step) with respect to the stress function O(r, t) by the method of inverse function. The time step was yielded by the Runge principle. The computational results in the form of frequency power spectra, wavelet spectra, and the surfaces of shell deflection for fixed time instances t = 300, t = 304, and t = 305 are reported in Table 3. The power spectrum obtained for n = 10 yields chaos spanned on five frequencies. In addition, the associated wavelet spectrum exhibits intermittency for n = 10, whereas for n = 20, chaos spanned on ω p occurs (For n = 40, one can detect ω p , ω p /2, ω p /4, and 3ω p /4.) Since the latter behavior is exhibited by both FFT and wavelet spectrum, where intermittency does not appear, the number of partitions of the radius plays a crucial role in identifying the scenario of transition from periodic to chaotic vibrations.
Furthermore, convergence of the deflection surfaces was investigated at the same time instances. As the reported graphs show, the surfaces coincide only for n = 120 and n = 140. To make the results more clear, each of the 3D surfaces is supplemented with the deflection isoclines. Figure 3a presents the time histories/signals in relation to the number n of partitions of the shell radius. (Figure 3b reports a window of the signal in the interval t ∈ [290; 330].) Remarkably, for n = 120 and n = 140, the graphs coincide.  Integration with respect to time was also carried out by both explicit and implicit Runge-Kutta methods of different orders, as it is shown in Table 1. It is worth mentioning that triangular or arbitrary matrices are used in explicit and implicit methods, respectively.
The methods RKF45, RKCK, and RK8pd allow for an automatic change of the integration steps as well as control of the integration errors. The carried out investigations showed that the fourth-order Runge-Kutta method is optimal to investigate the chaotic vibration of the considered flexible spherical shell subjected to harmonic load.
To discuss how the scenario to chaos will be changed while the number of nodes increases in the FDM, let us begin with the results reported in Table 2.
In this table the frequency power spectra are reported for the vibration signals obtained with the FDM for different node numbers n = 20; 40; 80; 120. For q 0 = 0.9, all frequency spectra exhibit ω p and ω p /2. Amplitudes of these harmonics are the same for all employed n. For q 0 = 0.1335, the regular bifurcations are exhibited on the frequencies ω p , ω p /2, and ω p /4 also for all n = 20; 40; 80; 120 nodes, but the values of harmonic amplitudes associated with ω p /4 are different for different n. This requires the use of higher number of nodes, i.e., the convergence of the results is achieved for n = 80; 120; 140; 160. In the case of q 0 = 0.1395, the vibration signals and their frequency characteristics strongly differ from each other: For n = 20, we have chaos with the distinguished ω p and ω p /2; for n = 40, we have also chaos, but with dominating ω p and ω p /4, whereas for n = 120; 140; 160, one deals with quasi-periodic vibration associated with the frequencies ω p , ω p /2, ω p /4, ω p /8, and ω p /16. For q 0 = 0.14, beginning already from n = 20, the obtained results coincide, exhibiting the system chaotic vibrations. Thus, the following conclusion can be yielded. A number of nodes employed in the FDM should be chosen differently for different vibrational regimes. Namely, periodic regimes require n = 20 nodes, whereas quasi-periodic/chaotic regimes  ω p require n = 80 or n = 120. It means that if one investigates a mechanical system with an infinite number of degrees of freedom, one has to investigate the convergence of numerical solutions obtained by the finite difference method in order to get reliable/true results. Table 3 contains frequency power spectra (FFT), 2D wavelet spectra, and deflection forms of the shell sur-faces at three different time instants, put one next to another for better clarity of the obtained results.  Table 4. It should be emphasized that to get reliable chaotic dynamics, a few different methods should be employed to compute LEs. In the present Table 3 Number of radius partitions n, power spectra, 2D wavelet spectra, and the shell deflection surfaces at three different time instants The carried out preliminary investigations allow one to perform reliable analysis of mechanical systems with infinitely many degrees of freedom.

Modes of vibrations (simple support)
In 1978, Feigenbaum [38] proposed a universal mechanism of period-doubling bifurcation route to chaos. The Feigenbaum scenario was validated for the simply supported (2.2) axially symmetric spherical shells subjected to transverse harmonic load in reference [39], where the following parameters were used: b = 4, ε = 0.01 (damping coefficient), ν = 0.3 (Poisson's ratio), t ∈ [0, 1400], ω p = ω 0 = 0.516, and q 0 ∈ [0.01, 0.15] (amplitude of harmonic excitation). As the excitation frequency ω p was equal to the shell frequency ω 0 , the case of the fundamental resonance was considered. Six Hopf bifurcations were detected and served later for estimation of the Feigenbaum constant yielding α 6 = 4.65608466. (The exact numerical value obtained by Feigenbaum equals 4.66916224 and differs by 0.28% from the above-mentioned numerical results.) The signals obtained in the present numerical experiments have the amplitude-time representation. However, in order to understand the behavior of any dynamic system under action of an arbitrary external load, and to control the system, one needs to have an insight into information hidden in the frequency space of the signal. For this reason, the wavelet transform approach was employed. Recall that the waveletbased analysis enables one to detect and trace the localized particular features of the signal, whereas the Fourier-based approach (FFT) makes it possible to follow the signal properties on the whole investigated time interval (see the description given in Sect. 3).
An increase in the number of new frequencies in the spectrum of vibration of the system yields a decrease in the informative properties of the employed wavelets, since harmonics with larger power shadow those with less power. Aside from this, it should be noted that the wavelet spectrum significantly depends also on the employed wavelet. (Each wavelet exhibits its own properties and peculiarities.) Therefore, to obtain reliable results, it is recommended to use both FFT and wavelet spectra simultaneously. The so far presented remarks concerning investigation of continuous structural systems were taken into account while constructing Table 5, where the Fourier frequency power spectra, 2D wavelet spectra, (w,ẇ,ẅ) phase plots, time history with marked time instants a, b, c, and the snapshots of the shell surface cross sections and the shell 3D shapes at a, b, c are reported for different values of q 0 . (A similar table was constructed when studying the problem with boundary conditions (2.3).) In the considered case, six period-doubling bifurcations were numerically detected and the Feigenbaum constant was estimated. The investigations showed that for the implemented boundary condition, a transition from periodic to chaotic vibrations is carried out following the Feigenbaum scenario. Remarkably, owing to the employed boundary conditions, one of the ribs at the surface of the shell deflections occurred regardless of the amplitude of the excitation load, i.e., the occurrence of the mentioned rib was not influenced by the number of the Hopf bifurcations. Furthermore, the mentioned rib appeared at the time instants when the shell top deflection was close to zero.

Modes of vibrations (rigid clamping)
In this section, vibrations of rigidly supported spherical shell (b = 8; boundary conditions (2.3)) subjected to the harmonic load q 0 = {0.01; 0.17; 0.21} and ω p = 1.34 close to ω 0 are studied. Namely, the case of the fundamental resonance is investigated. The associated dynamic characteristics, similar to those reported in Table 5, are shown in Table 6. In this case, a transition from regular to chaotic vibrations follows the Ruelle-Takens-Newhouse route, i.e., the increase in the excitation amplitudes yields the birth of the independent frequency, which is linearly dependent on the excita- Table 5 Fourier frequency power spectra, 2D wavelet spectra, (w,ẇ,ẅ) phase plots, time history with marked time instants a, b, c, shell surface cross sections, and shell 3D shapes at a, b, c for different values of q 0 (simple support)   tion frequency. For q 0 = 0.01, the shell exhibits a periodic regime. (Phase portraits present ellipses, whereas the LLE is less than zero.) The shell surface deflections exhibit two ribs, which occur at time instants when the shell center deflection is close to zero. An increase in the excitation load amplitude up to q 0 = 0.17 yields the occurrence of an independent frequency. Furthermore, the ribs are exhibited by the shell surface, and the LLE>0. The wavelet spectrum and the FFT spectrum also report the independent frequency, and a linear combination of the excitation and independent frequencies is observed. The investigations imply that in the case of the rigidly clamped axially symmetric spherical shell, a transition from regular to chaotic vibrations is realized via the Ruelle-Takens-Newhouse scenario along the shell contour. Moreover, deflection of the shell surfaces exhibits from one up to four ribs. It is worth mentioning that the Ruelle-Takens-Newhouse scenario plays a crucial role in the ribs occurrence. Namely, the presence of only one rib is associated with regular shell vibration, two ribs are associated with the occurrence of the independent frequency, whereas the occurrence of chaos yields the existence of up to three ribs.

Investigation of occurrence of ribs (simple nonmovable shell support)
This section is devoted to vibrations of the axially symmetric shell with simple nonmovable support of its contours and the shallow parameter b = 4. If the shell is subjected to transverse and uniformly disturbed harmonic load, a period-doubling route to chaos following the Feigenbaum scenario is detected and seven Hopf bifurcations can be estimated numerically. Contrary to the previous examples, the number of the detected ribs is not fixed, but changes in time (Table 7). Changes in the vibration character do not influence the number of the detected ribs, i.e., the Feigenbaum scenario does not influence the number of ribs. Table 6 Fourier frequency power spectra, 2D wavelet spectra, (w,ẇ,ẅ) phase plots, time history with marked time instants a, b, c, shell surface cross sections, and the shell 3D shapes at a, b, c for different values of q 0 (rigid clamping) ,   In what follows, vibrations of a spherical axially symmetric shell are considered for fixed parameters b = 5, ε = 0.1. The presence of the Sharkovskii's orbits is detected for the transverse uniform harmonic load q = q 0 sinω 0 t for the case of movable clamping (2.4) and for zero initial conditions (2.6). The shell is made from an elastic, homogenous, and isotropic material with ν = 0.3, and it is excited by the frequency ω p being equal to the fundamental shell frequency ω 0 = 0.773 (linear vibration). Phase portraits (3D), frequency power spectra, wavelet-based analysis results, and the shapes of the shell surface deflections are reported. Table 8 shows the Sharkovskii's ordering 3, 3 2 . The detected Sharkovskii's ordering occurred successively, and this novel result, to the knowledge of the authors, is reported for the first time in the literature. It should be mentioned that Sharkovskii's windows have been already detected and illustrated when studying vibration of beams.
For the applied boundary conditions, one timedependent rib was detected in the snapshots of the shell main cross sections. Here, no relationship between the vibration regimes and the number of ribs was found.

Concluding remarks
In the present paper, novel results of the study of chaotic vibration of flexible circular axially symmetric shallow shells subjected to sinusoidal transverse load were presented for four different boundary conditions . The study was conducted with the use of the finite difference method (FDM), which differentiates the study from the majority of research conducted with the finite element method (FEM).
The rigorous discussion of the advantages of the FDM can be found in numerous papers, including [40,41]. The shells considered in the present study Table 7 Fourier frequency power spectra, 2D wavelet spectra, (w,ẇ,ẅ) phase plots, time history with marked time instants a, b, c, and the snapshots of the shell surface cross sections and the shell 3D shapes at a, b, c for different values of q 0 (simple nonmovable shell support)   are axially symmetric, and hence, they are simpler for investigation in comparison with, for instance, geometric nonlinear shallow rectangular shells of the secondorder accuracy.
It should be emphasized that in the paper [42], the authors studied the dynamics of geometrically nonlinear beams by means of employing both the FEM in the Galerkin form and the FDM. The solutions obtained with the use of these methods coincided, but the FDM allowed one to obtain the solution 1.5 times faster than the FEM.
Finally, there are numerous products aimed at solving the application-oriented problems of the theory of elasticity, such as ANSYS, ABACUS, COMSOL, where FEM is used. One of the serious drawbacks of the mentioned commercial products is that they employ the dimensional equations and the computational time is very long.
The conducted investigations allow one to formulate the following general conclusions.

A general methodology for detection and analysis
of reliable/true chaos in the case on nonlinear vibration of flexible and axially symmetric shells sub-jected to uniformly disturbed harmonic load was proposed. 2. It was detected, illustrated, and discussed that the employed boundary conditions significantly influence the obtained solutions. For simple nonmovable and movable supports of the shell, for ω p = ω 0 , a transition from regular to chaotic shell vibration was realized following the Feigenbaum scenario, whereas for a rigidly clamped shell, a transition to chaos followed the Ruelle-Takens-Newhouse scenario. Finally, in the case of the movable clamped shell, a series of Sharkovskii's windows was reported. 3. It was found that the number of ribs exhibited by the shell surface strongly depends on the boundary conditions. The maximum number of ribs was detected for the shell with movable simple support, one rib was found for the shell with nonmovable simple support and the shell with movable clamping, whereas three ribs were detected in the case of rigid clamping. Table 8 Fourier frequency power spectra, 2D wavelet spectra, (w,ẇ,ẅ) phase plots, time history with marked time instants a, b, c, and the snapshots of the shell surface cross sections and the shell 3D shapes at a, b, c for different values of q 0 (movable clamping)