Application of the Lyapunov Exponents and Wavelets to Study and Control of Plates and Shells

1. Introduction

We study regular and chaotic vibrations of continuous mechanical systems using the following structural members: plates, flexible shells and cylindrical panels. They are often used in various measurement devices, and numerous ship and planes constructions. Owing to a high development of technology and due to industrial requirements, in many cases the mentioned structural members are subjected to action of high intensity loads being spatially and time dependent. The structural member exhibit either regular or chaotic dynamics, and hence an important question arises: How to predict safe and dangerous regimes of behavior of the studied mechanical objects?

In order to solve the mentioned problems, the investigations are carried out using the achievements of the qualitative analysis of differential equations and non-linear dynamics. Namely, we analyze time histories (signals), phase and modal portraits, Poincaré sections, autocorrelation functions, Lyapunov exponents, Fourier and wavelet spectra. Charts of vibration regimes versus the load excitation amplitude and frequency are constructed, which allow to control the vibration character of plates and shells.

Usually the data provided by numerical experiments are presented in time domain. In other words, we take time as an independent co-ordinate, and amplitude as a dependent co-ordinate, and the studied signal as analyzed through its amplitude-time representation. However, in order to understand deeply non-linear continuous systems subjected to various types of load actions and in order to fully understand the occurring dynamics, we have to apply the information hidden in the spectral signal characteristics. The Fourier transformation has been applied for a long time. However, it has been demonstrated recently that the Fourier analysis (FFT) is reliable only for the study of frequency components of stationary processes, i.e. the processes which through the whole period of investigation keep constant frequency components in time. It happens that in particular the dynamics of continuous mechanical systems may exhibit quite complicated output, and their frequency characteristics may change strongly in time. This is why in spite of the standard Fourier approach the wavelet analysis is applied allowing us to detect and understand many interesting non-linear phenomena of the mentioned mechanical systems.

Chaotic dynamics of structural members has been investigated by many researchers [1-10]. In this work we propose a novel approach to study non-linear vibrations of a plate based on the neural network approach and we analyze dynamics of flexible shells with constant stiffness and density subjected to harmonic load action. In the latter case mathematical model is built on the Kirchhoff-Love hypothesis and taking into account non-linear relation between deformation and displacement in the von Kármán form. This approach yields a system of non-linear PDEs regarding the deflection function and stresses (Airy’s function) as well as the system of equations regarding displacements [11]. We use further FDM with approximation $$O(h^2)$$ and BGM in higher approximations, which allows to study the system with infinite numbers of degrees of freedom without any truncation of the obtained system of ODEs, which is solved via the fourth-order Runge-Kutta method.

2. Lyapunov exponents computation via neural networks an the generalized Benettin’s algorithm

We illustrate and demonstrate the efficiency of the neural network approach to study a flexible plate with an infinite length. The governing equations are within the Kirchhoff hypothesis and they have the following non-dimensional form [11]:

$$\begin{array}{c}\frac{{\partial }^{2}u}{\partial x^2}+L_3(w,w)-\frac{{\partial }^{2}u}{\partial t^2}=0,\\ \frac{1}{{\lambda }^2}\left\{-\frac{1}{12}\frac{{\partial }^{4}w}{\partial x^4}+L_1(u,w)+L_2(w,w)\right\}+q-\frac{{\partial }^{2}w}{\partial t^2}-\epsilon \frac{\partial w}{\partial t}=0,\end{array}$$E1

where $L_1(u,w),L_2(w,w),L_3(w,w)$ – non-linear operators; $w(x,t)$– plate element bending in normal direction; $u(x,t)$– plate element longitudinal displacement; $\epsilon$– dissipation coefficient; $E$– Young modulus; h – height of the transversal panel cross section; $\gamma$– specific plate material gravity; g –Earth acceleration; t – time; $q=q_0\sin ({\omega }_{p}t)$– external load.

The non-dimensional parameters are as follows:

$$\begin{array}{l}\lambda =\frac{a}{h},\overline{w}=\frac{w}{h},\overline{u}=\frac{ua}{h^2},\overline{x}=\frac{x}{a},\overline{t}=\frac{t}{\tau },\\ \tau =\frac{a}{p},p=\sqrt{\frac{Eg}{\gamma }},\overline{\epsilon }=\frac{a}{p},\overline{q}=\frac{q{a}^4}{h^{4}E},\end{array}$$E2

and bars over the non-dimensional quantities have been already omitted in equations (1).

We demonstrate how to determine four first Lyapunov exponents applying pinned boundary conditions:

$w(0,t)=w(1,t)=u(0,t)=u(1,t)=w_{xx}^{//}(0,t)=w_{xx}^{//}(1,t)=0,$E3

and the following initial conditions

$$w(x,0)=\dot{w}(x,0)=u(x,0)=\dot{u}(x,0).$$E4

The boundary value problem (1), (3), (4) is reduced to the Cauchy problem via FDM (Finite Difference Method) of the second order accuracy. The obtained ODEs are solved by the Runge-Kutta method of the fourth and sixth orders. Validity and reliability of the obtained numerical results are confirmed by the FEM results (Finite Element Method). The initial problem of infinite dimension is substituted by that of finite dimension via partition of the interval $x\in [0,1]$ into 120 parts.

One of the ways to compute the spectrum of Lyapunov exponents is the neural network approach based on the generalized Benettin algorithm. It includes the following successive steps: 1. Choice of the appropriate time delay via tests; 2. Computation of an embedding space dimension; 3. Reconstruction of pseudo-phase trajectories using the method of time delays; 4. Neural network approximation; 5. Teaching of neural networks to compute successive iteration vectors; 6. Computation of the spectrum through the generalized Benettin algorithm with the help of the neural networks approach.

We apply the neural network with the following properties: It is an analog network regarding the input data (information is delivered through real numbers); It is self-organized with respect to its teaching aspects (output space of solutions is defined only through the input data); It belongs to the neural networks of straight signal distributions (all neural network couplings come from the input neurons and go to the output neurons); the neural network possesses dynamic couplings (control and improvement of synaptic couplings is carried out during the neural network self-teaching process $\left(\frac{dW}{dt}\ne 0\right),$ where W stands for the net weight coefficients). Let point $x_0$ belong to attractor A of a dynamical system. The trajectory of evolution of point is called the unperturbed trajectory. We choose the positive quantity $\epsilon$ essentially less than the attractor dimension. Next, we choose an arbitrary (perturbed) point in a way to satisfy the relation Then, we monitor the evolution of chosen and in time interval T, and the corresponding new points obtained after that time interval are denoted as and respectively. Vector is called the perturbation vector. We are ready to estimate λ:

${\tilde{\lambda }}_1=\frac{1}{T}\ln \frac{\Vert \Delta x_1\Vert }{\epsilon }.$E5

$q_0=0.125\times {10}^3$
a	b	c
d	e	f
$q_0=5\times {10}^3$
a	b	c
d	e	f
$q_0=7\times {10}^3$
a	b	c
d	e	f

Table 1.

Plate output characteristics

Time interval T is chosen in a way to keep the perturbation amplitude less than the linear dimensions of the space non-homogeneity as well as less than the attractor dimension. We consider the unit normalized perturbation vector and the corresponding new perturbation point ${\tilde{x}}_1^{/}=x_1+\Delta x_1^{/}.$ We extend the so far described approach by using points and ${\tilde{x}}_1$ instead of $x_0$ and respectively. We repeat the described procedure M times, and we may estimate λ as an average arithmetic quantity ${\tilde{\lambda }}_i$ of those obtained on each computation step. The proposed approach has been tested using the standard classical examples including that of the Henon map [12], the Lorenz system [13] and the logistic map.

We consider vibrations of our mechanical object with the following fixed parameters: $\lambda =50,$$\epsilon =1,$ ${\omega }_p=7,$ $q=q_0\sin ({\omega }_{p}t),$ and for the following amplitudes of the harmonic excitation: $q_0=0,125\cdot {10}^3;$$5\cdot {10}^3;$ $7\cdot {10}^3.$ In order to study chaotic dynamics of flexible plates we need to monitor and analyze the following output characteristics: time histories (a), phase (c) and modal portraits; phase portraits yielded by the neural networks approach (d); Fourier power spectra (b); wavelet spectra, Poincaré sections (e); spectra of Lyapunov exponents, where $$d$$ stands for the fractional part of dimension and $$h$$is the Kolmogorov-Sinai entropy (f); autocorrelation functions (some of them are reported in Table 1). Analysis of the obtained results implies that for $q_0=0,125\cdot {10}^3;$ periodic vibrations appear, whereas for $5\cdot {10}^3$ chaos is exhibited, and for $7\cdot {10}^3$ the hyper-chaotic vibrations occur.

3. Wavelets approach to study plate dynamics

One of the first important tasks to be solved is associated with the choice of a suitable wavelet, which fits well with the stated problem. In order to solve this query we have studied the non-stationary signal (Table 2) obtained via the numerical experiment as an output of the mathematical model of the rectangular flexible isotropic plate subjected to the periodic shear load acting in the shell volume unit. The mathematical model is as follows [11]:

$\begin{array}{l}\frac{1}{12(1-{\mu }^2)}\left({\nabla }_{\lambda }^{4}w\right)-L(w,F)+\frac{{\partial }^{2}w}{\partial t^2}+\epsilon \frac{\partial w}{\partial t}-q(x_1,x_2,t)+2S\frac{{\partial }^{2}w}{\partial x_1\partial x_2}=0,\\ {\nabla }_{\lambda }^{4}F+\frac{1}{2}L(w,w)=0,\end{array}$E6

where: $${\nabla }_{\lambda }{}^4=\frac{1}{{\lambda }^2}\frac{{\partial }^4}{\partial x_1{}^4}+{\lambda }^2\frac{{\partial }^4}{\partial x_2{}^4}+2\frac{{\partial }^4}{\partial x_1{}^2\partial x_2{}^2}$$, $$L(w,F)=\frac{{\partial }^{2}w}{\partial x_1{}^2}\frac{{\partial }^{2}F}{\partial x_2{}^2}+\frac{{\partial }^{2}w}{\partial x_2{}^2}\frac{{\partial }^{2}F}{\partial x_1{}^2}-2\frac{{\partial }^{2}w}{\partial x_1\partial x_2}\frac{{\partial }^{2}F}{\partial x_1\partial x_2}$$– the known non-linear operator, whereas w and F stand for the plate deflection and Airy’s function, respectively.

System (3.1) is reduced to the non-dimensional form using the following non-dimensional parameters: $\lambda =a/b,$$x_1=a{\overline{x}}_1,$ $x_2=b{\overline{x}}_2$ – non-dimensional parameters regarding $$x_1$$ and $$x_2$$, respectively; $$w=2h\overline{w}$$– deflection; $F=E{(2h)}^3\overline{F}$– Airy’s function; $t=t_0\overline{t}$– time; $q=\frac{E{(2h)}^4}{a^2{b}^2}\overline{q}$– external load; $\epsilon =(2h)\overline{\epsilon }$– damping coefficient of the surrounding medium, $S=\frac{E{(2h)}^3}{ab}\overline{S}$– external shear load. In the equations bars have already been omitted over the non-dimensional quantities. The following notation is introduced: $$a,b$$– plate dimensions regarding $$x_1$$ and $$x_2$$, respectively; $\mu$– Poisson’s coefficient. Zero value initial conditions and the following boundary value conditions are attached to system (2.6):

$\begin{array}{l}w=0,\frac{{\partial }^{2}w}{\partial x_1^2}=0,F=0,\frac{{\partial }^{2}F}{\partial x_1^2}=0for{x}_1=0;1;\\ w=0,\frac{{\partial }^{2}w}{\partial x_2^2}=0,F=0,\frac{{\partial }^{2}F}{\partial x_2^2}=0for{x}_2=0;1.\end{array}$E7

The external harmonic shear load has the form $S=s_0\sin {\omega }_{p}t.$ PDEs governing dynamics of our investigated plate are reduced to the ODEs via the FDM (Finite Difference Method) with the approximation О(h²) regarding spatial co-ordinates. Next, ODEs are solved via the fourth-order Runge-Kutta method, and additionally on each of the iterations a large system of linear algebraic equations should be solved with respect to the stress (Airy’s) function. Time integration step has been chosen using the Runge rule. The partition number of spatial co-ordinates is n=14 while applying FDM. Validity and reliability of the obtained results regarding the number of partitions has been discussed by Awrejcewicz et al. [14].

The term $2S\frac{{\partial }^{2}w}{\partial x_1\partial x_2}$ introduced in the governing equations exhibits the action of shear stresses located in shell middle plane and it essentially influences non-linear dynamics of the investigated shell. The numerical simulation indicates that the output signal (shell vibrations) may change in time repeatedly. We apply this signal to choose a methodology suitable for the investigation of non-stationary processes, and in addition, we illustrate advantages and disadvantages of the standard Fourier approach versus the wavelet transform procedure. The studied signal has been obtained analyzing the system with the following fixed parameters: $s_0=8.4$ and ${\omega }_p=26.$ We show that frequency characteristics taken in different time intervals essentially differ from each other. It should be emphasized that the system stability loss occurs not only via the change of chosen control parameters but also keeps all of them fixed owing to the system time evolution. In the first time interval $t\in$ [50;56] the shell exhibits two frequency quasi-periodic vibrations. Instead of the vanished excitation frequency, two independent frequencies have appeared. A further long time evolution of chaotic vibrations with the exhibition of a few dependent frequencies is observed. In the Fourier spectrum the excitation frequency is not visible. The last studied time interval corresponds to harmonic vibrations, which is also in agreement with the Fourier power spectrum.

We constructed wavelet spectra regarding the mentioned signal. We applied the following wavelets: Haar, Shannon-Kotelnikov, Meyer, Daubechies wavelets from db2 up to db 16, Coiflets and symlet wavelets, as well as the Morlet and Gauss (real and imaginary) wavelets, on the basis of the derivatives from 2 to 16. Haar and Shannon-Kotelnikov wavelets are not suitable for the analysis of shell structures. The first one is badly localized regarding frequency, whereas the second one, contrary to the previous wavelets, is badly localized in time. On the other hand, the analysis regarding the Doubechies wavelets, as well as symlet and Coiflets wavelets implies an increase of the frequency resolution assuming that the filter properties are increased. Neglecting differences regarding the wavelet forms and the associated filters, the wavelet spectra obtained through the Dobechies wavelets as well as symlet and Coiflets wavelets are practically identical. However, their localization with respect to frequency is not suitable for the analysis of continuous systems dynamics. In the case of the Gauss functions, an increase of their derivative order implies an increase of the frequency resolution.

Fourier spectrum (50<t<56) of quasi-periodic vibrations	Fourier spectrum (130<t<225) of chaotic vibrations	Fourier spectrum (230<t<286) of periodic vibrations
Signal

Table 2.

Fourier spectra and a signal

Table 3 gives results associated with the application of different wavelets (Meyer, Morlet, complex Morlet, real and complex Gauss with 16 derivative order, Daubechies) to analyze non-linear shell vibrations. One may conclude from Table 3 that the localization regarding frequency increases with an increase of the number of the wavelet zero moments.

4. Equations governing dynamics of flexible elastic shells

In a frame of non-linear classical theory we study a shell with constant stiffness and density and harmonic load $$q=q_0\sin ({\omega }_{p}t)$$, where $$q_0$$ - amplitude of excitation, $${\omega }_p$$- frequency of excitation. In the rectangular co-ordinates the 3D space is:

$$\Omega =\{x_1,x_2,x_3|(x_1,x_2)\in [0;a]\times [0;b],x_3\in [-h/2;h/2]\},$$$$0\le t<\infty$$.

Wavelet spectrum: Меyеr	Wavelet spectrum: Daubechies 16	Wavelet spectrum: complex Morlet
Wavelet spectrum: Morlet	Wavelet spectrum: Gauss 16	Wavelet spectrum Gauss: 16 complex

Table 3.

Wavelets spectra

The governing non-dimensional equations are given in the hybrid form:

$$\begin{array}{l}\frac{1}{12(1-{\mu }^2)}\left({\nabla }_{\lambda }{}^{4}w\right)-k_{x_2}\frac{{\partial }^{2}F}{\partial x_1{}^2}-k_{x_1}\frac{{\partial }^{2}F}{\partial x_2{}^2}-L(w,F)+\frac{{\partial }^{2}w}{\partial t^2}+\epsilon \frac{\partial w}{\partial t}-q(x_1,x_2,t)=0,\\ {\nabla }_{\lambda }{}^{4}F+k_{x_2}\frac{{\partial }^{2}w}{\partial x_1{}^2}+k_{x_1}\frac{{\partial }^{2}w}{\partial x_2{}^2}+\frac{1}{2}L(w,w)=0,\end{array}$$E8

whereas equations regarding displacements follow

$$\begin{array}{l}\frac{{\partial }^{2}u}{\partial x_1^2}+\frac{1-\mu }{2}\frac{{\partial }^{2}u}{\partial x_2^2}+\frac{1+\mu }{2}\frac{{\partial }^{2}v}{\partial x_1\partial x_2}+\frac{\partial w}{\partial x_1}\frac{{\partial }^{2}w}{\partial x_1^2}-\left(k_{x_1}+\mu k_{x_2}\right)\frac{\partial w}{\partial x_1}+\frac{1+\mu }{2}\frac{\partial w}{\partial x_2}\frac{{\partial }^{2}w}{\partial x_1\partial x_2}+\frac{1-\mu }{2}\frac{\partial w}{\partial x_1}\frac{{\partial }^{2}w}{\partial x_1^2}-\frac{{\partial }^{2}u}{\partial t^2}=0,\\ \frac{1+\mu }{2}\frac{{\partial }^{2}u}{\partial x_1\partial x_2}+\frac{{\partial }^{2}v}{\partial x_2^2}+\frac{1-\mu }{2}\frac{{\partial }^{2}v}{\partial x_1^2}-\left(\mu k_{x_1}+k_{x_2}\right)\frac{\partial w}{\partial x_2}+\frac{\partial w}{\partial x_2}\frac{{\partial }^{2}w}{\partial x_2^2}+\frac{1+\mu }{2}\frac{\partial w}{\partial x_1}\frac{{\partial }^{2}w}{\partial x_1\partial x_2}+\frac{1-\mu }{2}\frac{\partial w}{\partial x_2}\frac{{\partial }^{2}w}{\partial x_1^2}+\frac{{\partial }^{2}v}{\partial t^2}=0,\\ \frac{1}{{\lambda }^2}{[{\nabla }^{4}w-\left(k_x+\mu k_{x_2}\right)\frac{\partial u}{\partial x_1}-\left(\mu k_{x_1}+k_{x_2}\right)\frac{\partial v}{\partial x_2}+\left(k_{x_1}^2+k_y^2+2\mu k_{x_1}{k}_{x_2}\right)w-\frac{k_{x_1}+\mu k_{x_2}}{2}\left(\frac{\partial w}{\partial x_1}\right)}^2-\\ -\frac{k_{x_2}+\mu k_{x_1}}{2}{\left(\frac{\partial w}{\partial x_2}\right)}^2-\frac{\partial }{\partial x_1}\left\{\frac{\partial w}{\partial x_1}\left[\frac{\partial u}{\partial x_1}+\mu \frac{\partial v}{\partial x_2}\right]+\frac{1-\mu }{2}\frac{\partial w}{\partial x_2}\left(\frac{\partial u}{\partial x_2}+\frac{\partial v}{\partial x_1}\right)\right\}-\\ -\frac{\partial }{\partial x_2}\left\{\frac{\partial w}{\partial x_2}\left[\frac{\partial v}{\partial x_2}+\mu \frac{\partial u}{\partial x_1}\right]+\frac{1-\mu }{2}\frac{\partial w}{\partial x_1}\left(\frac{\partial u}{\partial x_2}+\frac{\partial v}{\partial x_1}\right)\right\}]-q+{\epsilon }_1\frac{\partial w}{\partial t}+\frac{{\partial }^{2}w}{\partial t^2}=0,\end{array}$$E9

where: $${\nabla }_{\lambda }{}^4=\frac{1}{{\lambda }^2}\frac{{\partial }^4}{\partial x_1{}^4}+{\lambda }^2\frac{{\partial }^4}{\partial x_2{}^4}+2\frac{{\partial }^4}{\partial x_1{}^2\partial x_2{}^2}$$, $$L(w,F)=\frac{{\partial }^{2}w}{\partial x_1{}^2}\frac{{\partial }^{2}F}{\partial x_2{}^2}+\frac{{\partial }^{2}w}{\partial x_2{}^2}\frac{{\partial }^{2}F}{\partial x_1{}^2}-2\frac{{\partial }^{2}w}{\partial x_1\partial x_2}\frac{{\partial }^{2}F}{\partial x_1\partial x_2}$$– known non-linear operator, $$w(x_1,x_2,t)$$– element normal; $$u(x_1,x_2,t)$$, $$v(x_1,x_2,t)$$– element displacement regarding $$x_1$$ and $$x_2$$, respectively; $$F(x_1,x_2,t)$$- stress function; $$\epsilon$$ – dissipation coefficient; $$E$$– Young modulus; $$h$$ – height of the transversal panel cross section; $$\gamma$$– specific unit gravity of the shell material; $$g$$– Earth acceleration; $$t$$– time; $$q=q_0\sin (\omega t)$$ – external load.

The following non-dimensional parameters are introduced:

$$\begin{array}{l}\lambda =\frac{a}{h},{\lambda }_1=\frac{a}{b},\overline{w}=\frac{w}{h},\overline{u}=\frac{ua}{h^2},\overline{v}=\frac{va}{h^2},\overline{x_1}=\frac{x_1}{a},\overline{x_2}=\frac{x_2}{b},\overline{t}=\frac{t}{\tau },\tau =\frac{a}{c},\\ c=\sqrt{\frac{Eg}{\gamma }},\overline{\epsilon }=\frac{a}{c},\overline{q}=\frac{q{a}^4}{h^{4}E},\overline{k_{x_1}}=\frac{k_{x_1}a}{\lambda },\overline{k_{x_2}}=\frac{k_{x_2}b}{\lambda },\overline{F}=\frac{F}{E{h}^3}.\end{array}$$E10

System of differential equations (4.1)-(4.2) should be supplemented by boundary and initial conditions (see [11]. Since we cannot solve the stated problems analytically, we reduce the problem to ODEs and solve it numerically by the fourth-order Runge-Kutta method (see [15] for more details).

As the first example, we consider a spherical rectangular shell governed by equations (4.1) and with the following homogeneous conditions:

$$w{|}_{\Gamma }=0;M_n{|}_{\Gamma }=0;N_n{|}_{\Gamma }=0;{\epsilon }_n{|}_{\Gamma }=0\text{for}{x}_1=0;1,x_2=0;1,$$E11

which can be recast to the form

$$w=0;\frac{{\partial }^{2}w}{\partial x_1{}^2}=0;F=0;\frac{{\partial }^{2}F}{\partial x_1{}^2}=0;w=0;\frac{{\partial }^{2}w}{\partial x_2{}^2}=0;F=0;\frac{{\partial }^{2}F}{\partial x_2{}^2}=0$$E12

and the following initial conditions

$$w(x_1,x_2){|}_{t=0}=0,\begin{array}{c}\end{array}\frac{\partial w}{\partial t}=0.$$E13

Geometric parameters of the shell curvature $$k_{x_1}=k_{x_2}=24$$, $$\lambda =1$$, and the damping coefficient $$\epsilon =1$$. We apply the BGM in Vlasov’s form, and the we are looping for the functions w and F, satisfying (4.5), in the following form

$$w={\displaystyle \sum _{i=1}^N{\displaystyle \sum _{j=1}^N{A}_{ij}(t)\sin (i\pi x_1)\sin (j\pi x_2),}}F={\displaystyle \sum _{i=1}^N{\displaystyle \sum _{j=1}^N{B}_{ij}(t)\sin (i\pi x_1)\sin (j\pi x_2)}}.$$E14

5. Numerical analysis of shells non-linear dynamics

In what follows we investigate vibrations of flexible shells via the BGM versus the partition number $$N$$ in periodic (Fig. 1) and chaotic (Fig. 2) zones. We consider the point $$A(q_0,{\omega }_p)=A(5;25)\in \left\{q_0,{\omega }_p\right\}$$ marked on the type vibration chart (Fig. 3), which belongs to harmonic vibrations zone. For all $$N$$ harmonic (one frequency) vibrations are observed. Increase of the approximations number $$N=11,13,15$$ implies the full coincidence regarding both frequency and amplitude. Further, we study the BGM convergence in a chaotic zone (Fig. 2). We analyze the point $$B(q_0,{\omega }_p)=B(180;38)\in \{q_0,{\omega }_p\}$$ on the chart (Fig. 3). Although we cannot achieve the signal convergence as it occurred in the previous case, but we get the convergence in average sense through Fourier integrals of power spectra.

Figure 1.

Time histories $$(w(0.5,0.5,t))$$ and power spectrum for different approximations of the Bubnov-Galerkin method for $$t\in [40;40.4]$$.

Figure 2.

Time histories $$(w(0.5,0.5,t))$$ and power spectrum for different approximations of the Bubnov-Galerkin method for $$t\in [40;41]$$.

Figure 3.

Vibration charts in the control parameters plane $$\{q_0,{\omega }_p\}$$.

We investigate further the convergence of the FDM versus a number of partitions of the mesh $$m\times n$$ for a flexible rectangular shell taking into account the same parameters. We consider first the shell center motion in an harmonic zone (Fig. 4). Let us fix the point $$A(q_0,{\omega }_p)=A(5;25)\in \left\{q_0,{\omega }_p\right\}$$ on the vibration chart. Three curves are reported in Fig. 4 for $$n=m=8;16;32$$. For all partition points $$n=m=8;16;32$$ harmonic vibrations are observed. On the other hand for given $$n=m$$ time histories differ from each other (see Fig. 5), but their convergence regarding the Fourier power spectra are evident ($$n=m=16;32$$).

Figure 4.

Deflections $$(w(0.5,0.5,t))$$ and power spectra for different partitions $$m=n=8;16;32$$ in a periodic zone for $$t\in [40;40.5]$$.

Figure 5.

Dependence $$(w(0.5,0.5,t))$$ and power spectra for different partition $$m=n=8;16;32$$ in a chaotic zone for $$t\in [40;45]$$.

We now compare the results obtained through two qualitatively different approaches, i.e. FDM and BGM. The convergence of those two methods is numerically confirmed with respect to time histories and Fourier power spectra for small amplitude of excitation. Although in the case of chaotic vibrations the convergence regarding time series is not achieved, but it is achieved with respect to integral Fourier characteristics.

In what follows we address the problem of vibration chart identification versus a step of variation of the control parameters $$q_0$$ and $${\omega }_p$$. We consider the vibration charts for five applied resolutions given in Table 4: а – resolution $$50\times 50$$, b -$$100\times 100$$, c –$$200\times 200$$, d –$$300\times 300$$, e -$$400\times 400$$.

$$50\times 50$$	$$100\times 100$$	$$200\times 200$$
$$300\times 300$$	$$400\times 400$$

Table 4.

Charts of shell vibrations in the $$\{q_0,{\omega }_p\}$$ plane.

Increase of the resolution implies the results improvement. Results obtained for the cases 1d and 1e coincide in full. In Figure 3 vibration charts for the rectangular spherical shell, which have been obtained using the BGM for N = 11, as well as the FDM for partition number 8х8 and accuracy $$O(h^2)$$ regarding the spatial coordinate are shown. In both cases the initial problem is solved via the fourth order Runge-Kutta method. It should be emphasized that in order to remove potential errors in the obtained results and in order to confirm its reliability and validity, one needs to apply different numerical approaches while studying non-linear dynamics of 2D mechanical objects. The obtained vibration charts allow controlling the vibration regimes and a transition form dangerous zones to those of the required/engineering acceptable zones. One may observe that both charts obtained via the qualitatively different approaches are close to each other. For small values of the load ($$q_0\le 30$$) the system exhibits periodic vibrations. Increase of the excitation amplitude implies occurrence of chaotic dynamics. There are also zones of the Hopf bifurcations and two-frequency quasi-periodic orbits. However, both charts exhibit a qualitative difference in the resonance zone ($${\omega }_p={\omega }_0$$), i.e. errors introduced by numerical techniques increase in the resonance conditions. However, increase of of the partition number of the applied mesh as well as increase of the applied terms of the series for the case of BGM again yields the reliable and validated results, although with a higher costs of time computation.

Approximation of $$N=15$$ for the BGM and of $$n=m=16$$for the FDM are most suitable to keep the validated results as well as economically reasonable computation time. Computational time of the BGM is lower than that of the FDM, since in the latter case we need to solve a system of algebraic equations (for $$n=m=8$$ we have 64 equations, for $$n=m=16$$-256 equations, $$n=m=32$$-1024 equations), which requires additional computational time.

Therefore, in the case of our 2D mechanical object we have achieved only the integral convergence regarding the Fourier spectrum, and hence one my study a 1D problem.

The second computational example deals with the infinite cylindrical panel harmonically and transversally loaded. In this case equations (4.2) yield

$$\begin{array}{c}\frac{{\partial }^{2}u}{\partial x_1{}^2}-k_{x_1}\frac{\partial w}{\partial x_1}+L_3(w,w)-\frac{\partial u^2}{\partial t^2}=0,\\ \frac{1}{{\lambda }^2}\left\{-\frac{1}{12}\frac{{\partial }^{4}w}{\partial x_1{}^4}+k_{x_1}\left[\frac{\partial u}{\partial x_1}-k_{x_1}w-\frac{1}{2}{\left(\frac{\partial w}{\partial x_1}\right)}^2-w\frac{{\partial }^{2}w}{\partial x_1{}^2}\right]+L_1(u,w)+L_2(w,w)\right\}+q-\frac{{\partial }^{2}w}{\partial t^2}-\epsilon \frac{\partial w}{\partial t}=0.\end{array}$$E15

with the following boundary

$$w(0,t)=w(1,t)=u(0,t)=u(1,t)=w\text{'}{\text{'}}_{x_1{x}_1}(0,t)=w\text{'}{\text{'}}_{x_1{}_{}{x}_1}(1,t)=0,$$E16

and initial conditions

$$w(x_1,0)=\dot{w}(x_1,0)=u(x_1,0)=\dot{u}(x_1,0)=0.$$E17

	$50\times 50$	$100\times 100$	$200\times 200$	$300\times 300$
$n=16$
$n=32$
$n=40$

Table 5.

Vibration charts $$\{q_0,{\omega }_p\}$$ for different n and different resolutions

n=16	n=32	n=60	n=80

Table 6.

Time histories $$(w(0.5,t))$$, Fourier spectra, and Poincaré maps for different panel partitions.

Observe that neither the input governing equations (5.1) nor more complex system of equations (4.2) cannot be solved analytically. The results obtained through FDM have been compared with the results obtained FEM (Finite Element Method) in the Bubnov-Galerkin form. The following characteristics are applied: time history $$(w(0.5,t))$$, power spectrum and the Morlet wavelets. Charts of the vibration regimes present all peculiarities and a general picture of the studied non-linear process within the investigated intervals of the control parameters. However, one would expect an optimal choice of the applied nodes ($$n$$) of the mesh, as well a number of partition of the excitation frequency ($${\omega }_p$$) and the excitation amplitude ($$q_0$$). In Table 5 charts of the vibration regimes on the control parameters plane $$\{{\omega }_p,q_0\}$$ and the interval partition regarding $$x_1\in [0;1]$$, including the applied color notation is presented for $$\{{\omega }_p,q_0\}=50\times 50;100\times 100;200\times 200;300\times 300$$, and the interval $$x_1\in [0;1]$$ is divided into 16, 32 and 40 parts.

We have already mentioned that a crucial role in analysis plays the time of computation to get time histories and carry out their analysis. Taking into account the reported results we conclude that the most optimal chart resolution is that of $$200\times 200;$$ $$300\times 300.$$

It should be noted that the constructed charts on a basis of only power spectra are not sufficient to decide about the results convergence versus the partition number ($$n$$) regarding the spatial coordinate. In order to get the validated results we apply the following fixed parameters (boundary conditions (5.2) and initial conditions (5.3)):

$$\epsilon =1;$$ $$\lambda =\frac{a}{h}=50;$$ $${\omega }_p=8,625;$$ $$q_0=59000;$$ $$k_{x_1}=0.$$

For the given parameters the system is in a chaotic regime, what is approved by the reported charts (Table 5, point А). In what follows we analyze the obtained signal versus number of partition of the spatial coordinate. Namely, we construct the power spectra and Poincaré maps (Table 6) for n=16, 32, 40.

Comparison of the amplitudes of the analyzed time histories allows to conclude that the convergence is achieved for. Now, taking into account the Poincaré maps we that for $$n\ge 60$$ a chaotic strange attractor appears, which has not been detected for a smaller number of applied partitions. The Fourier power spectrum exhibits a remarkable localization of dominating frequencies including the excitation frequency for different introduced partitions, but the convergence is not achieved.

Therefore, taking into account the mentioned remarks, the following parameters have been applied through computations:

• PDEs are solved via FEM;

• Charts resolution –$$200\times 200;$$

• Number of partitions regarding the spatial coordinate –$$n=80$$.

In the chapter part devoted to study non-linear vibrations of shells, we have chosen suitable charts resolutions and partition numbers to follow various scenarios of transitions form regular to chaotic dynamics and to optimal reduction of the computational time simultaneously keeping the reliable and validated results.

6. Concluding remarks

In the chapter part devoted to plate analysis we have shown that the obtained validated spectra of the Lyapunov exponents allow for the estimation of Kaplan-Yorke dimension, Sinai-Kolmogorov entropy, and velocity of the phase space compression. Furthermore, we have illustrated that the complex Morlet and Gauss wavelets have better localization with respect to frequency, in comparison to their real analogues, but time localization is better for the real wavelets. Therefore, one may either apply real or complex Morlet and Gauss (order bigger than 16) wavelets while studying plates/shells dynamics.

Acknowledgments

This work has been supported by the grant RFFI No 12-01-31204, and the National Science Centre of Poland under the grant MAESTRO 2, No. 2012/04/A/ST8/00738, for years 2013-2016.