Parametric vibrations of graphene sheets based on the double mode model and the nonlocal elasticity theory

Parametric vibrations of the single-layered graphene sheet (SLGS) are studied in the presented work. The equations of motion govern geometrically nonlinear oscillations. The appearance of small effects is analysed due to the application of the nonlocal elasticity theory. The approach is developed for rectangular simply supported small-scale plate and it employs the Bubnov–Galerkin method with a double mode model, which reduces the problem to investigation of the system of the second-order ordinary differential equations (ODEs). The dynamic behaviour of the micro/nanoplate with varying excitation parameter is analysed to determine the chaotic regimes. As well the influence of small-scale effects to change the nature of vibrations is studied. The bifurcation diagrams, phase plots, Poincaré sections and the largest Lyapunov exponent are constructed and analysed. It is established that the use of nonlocal equations in the dynamic analysis of graphene sheets leads to a significant alteration in the character of oscillations, including the appearance of chaotic attractors.


Introduction
It is clear that nanotechnology has taken a large place in modern industry in recent years. This is due to the fact that micro and nanoobjects have excellent mechanical, chemical, electrical, and magnetic properties. Such objects are widely used as elements of the micro(nano)electromechanical systems (MEMS and NEMS), resonators, sensors, atomic force microscopes, DNA detectors, thin films, and others [1][2][3][4][5][6][7]. However, the experimental and theoretical studies carried out have shown that the classical (sizeindependent) theory in the investigation of objects with dimensions in nano/microscale cannot guarantee correct outcomes. That is why size-dependent continuum theories have been developed and applied to various problems. Mention should be made of the micropolar elasticity theory by the Cosserat brothers [8], as well as the couple stress theory by Mindlin and Tiersten, Toupin, Koiter [9][10][11]. Applied in our paper the nonlocal elasticity theory by Eringen [12] uses the fact that the stress at a given point depend on strains at all other points in the body and introduces the supplementary nonlocal parameter. It should also be noted the strain gradient theory established by Lam et al. [13] including three material length scale parameters, and the modified couple stress theory by Yang et al. [14], which contains only one additional material length scale parameter, moreover, in this case there is a symmetric couple stress tensor.
As mentioned above, the proposed work uses the nonlocal theory of elasticity. There are many works in which linear vibrations and buckling of isotropic and orthotropic micro and nanoplates are studied [15][16][17][18][19]. A significantly smaller number of works is devoted to solving problems in a geometrically nonlinear formulation with various types of loading as well as under magnetic and electric excitation [20][21][22][23]. Moreover, in these works, periodic regimes are studied based on the single-mode approximation of the deflection. At the same time, numerous works based on the classical theory have shown that the use of two or more mode models makes it possible to detect chaotic regimes [24,25]. It is also worth mentioning here the works where the nonlocal effects taken into account. Nonlinear modes of vibration of single-layered graphene sheets are studied by Ribeiro and Chuaqui [26] by Bubnov-Galerkin approach as well as harmonic balance method. Awrejcewicz et al. [27] analyse the forced nonlinear vibrations of nano/micro plates based on multiple scale method Furthermore, when the plate is subjected to periodic in-plane load the parametric resonance may occur [28,29] and this matter should be investigated thoroughly.
A literature review showed that this problem of parametric vibrations of micro/nanoplates remains insufficiently studied. Those few recent studies of micro nanostructures use the classical Bolotin method for dynamical stability analysis. Gholami et al. studied axial buckling and dynamic stability of functionally graded microshells in [20]. The dynamic stability of Timoshenko nanobeams subjected to an axial load studied in the work of Behdad et al. [30], Bolotin's approach applied to the investigation of dynamic stability of functionally graded microbeams by Ke and Wang [31], where authors used the modified couple stress theory.
In our work we considered the orthotropic smallscale plate compressed by harmonic load in-plane along the axis Ox. The approach used in this work is based on the Bubnov-Galerkin method, while the solution is presented in the Navier form, taking into account two modes. This helped us to reduce the resolving size-dependent partial differential system to ODEs the coefficients of which contain a nonlocal parame- ter. Integration of the obtained system is carried out by the Runge-Kutta method. The complex behavior was investigated by varying the load parameter and nonlocal parameter in order to detect chaotic oscillations and analyse the transition to them.
The manuscript has five parts. The first section is aimed at the introduction, the next part contains the formulation of the problem based on the nonlocal elasticity theory. Following the statement part, the third section is devoted to describing the method used Bubnov-Galerkin approach and double mode model. The fourth part contains the results of numerical simulations and bifurcation analysis. The paper ends with some conclusions.

Mathematical formulation
The presented work is devoted to the investigation of the parametric vibrations of SLGS. The equations governing the problem use the nonlocal elasticity theory as well as are based on Von Kármán nonlinear theory and Kirchhoff-Love Hypothesis. The nanoplate is subjected to the in-plane uniaxial periodic force, see Fig. 1.
In contrast the classical theory, according to the nonlocal elasticity theory the constitutive relation for the nonlocal stress tensor at a point x is presented in the following form In formula (1) σ , σ are nonlocal and local stress tensors, K |X − X |, τ is the nonlocal modulus, τ = e 0 α/l, where α is the internal characteristic length (distance between C-C bonds, lattice parameter) and l stands for external characteristic length, whereas e 0 is a constant appreciating to material for adjusting to experimental results or molecular dynamics results [12,32].
The system of the movement equation is presented in the mixed form [33,34], for this, the stress function F is introduced in the following way where N x , N y , N xy are resultant in-plane forces. Note that the problem considered assumes neglecting of propagation of elastic waves in longitudinal equations. A detailed derivation of the governing system relating the stress function F and plate deflection w one can find in [22,23]: where , E 1 , E 2 are Young's moduli, ν 12 , ν 21 are Poisson's ratios, G 12 is shear modulus, h and ρ are thickness and density of the plate, whereas ε stands for the damping coefficient, μ is nonlocal parameter, introduced by the formula, also Function p(t) is periodic load acting in plane along Ox axis. The compatibility equation for strains [33,34] in the middle surface allows to obtain the next equation where The study was carried out for the case of a simply supported plate with movable edges:

Bubnov-Galerkin approach with double mode model
Let us present the deflection of the plate w(x, y, t) in the Navier form considering two vibration modes as follows: where w 1 , w 2 are bi-modal amplitudes, sin π x a sin π y b and sin 2π x a sin 2π y b are shape functions. The presentation of the stress function follows from the relation (6), where the expression for the deflection (9) is substituted in the right-hand side: thus, the stress function takes the form with coefficients defined by formulas , , Relations (8) allows to get coefficients p 1 = 0, p 2 = 0. Taking into account (9), (11) and applying the Bubnov-Galerkin method where yields the system of ordinary differential equations: with coefficients We consider the in-plane periodic loading applied along two opposite edges in the following form p x = p 0 cos ωt, where p 0 and ω are excitation amplitude and frequency. Using the dimensionless ratios one can get to the system of equations y 1 + δy 1 + ξ 0 y 1 − ξ p y 1 cos 2τ + ξ 1 y 1 y 2 2 + ξ 2 y 3 1 = 0, y 2 + δy 2 + η 0 y 2 − η p y 2 cos 2τ + η 1 y 2 y 2 1 + η 2 y 3 2 = 0.
The coefficients of the obtained system are defined by the formulas

Numerical results
Numerical results were obtained for a plate made of isotropic material with the following properties [15]: Note that for an isotropic material it follows It should be aware of in the literature there is a very limited number of works devoted to nonlinear vibrations of micro/nanoplates, so the comparison of the results was carried out for free linear vibrations with results published in [15], this comparative study can be viewed in our previous work [35]. It is also worth noting that in the case of an isotropic material within the framework of the classical theory (the nonlocal parameter μ = 0) and without taking into account the in-plane periodic load, we get system of nonlinear second-order ordinary differential equations where the obtained coefficients coincide with those presented in the work [24], which also confirms the results of our work. The numerical simulations are performed for mentioned material properties (20), they are assumed to be fixed as well as the damping parameter δ = 1. For studying parametrical resonant case we took the excitation frequency is equal to the double first natural frequency of plate ω = 2 √ α 0 . Thus, the dynamical behaviour of the system, chaotic vibrations appearance and the manifestation of small-scale effects are analysed for varying in-plane load parameter ξ p (η p = 4ξ p ) and the nonlocal parameter μ.
For numerical solution of ordinary differential equations we have used the package Wolfram Mathematica 12.0 and the built-in function NDSolve with default options. Therefore the integration methods and time steps were selected automatically. The bifurcation diagrams and Poincaré sections presented below were made on the basis of stroboscopic sampling of the system state at times satisfying the condition τ = πi (i ∈ Z ). Each bifurcation plot is composed of 600 Poincaré sections performed for 600 constant values of the bifurcation parameter, evenly distributed over the range of changes of this parameter. Each Poincaré map included in the bifurcation plot consists of 300 points corresponding to the steady state motion after omitting the initial 100 points of the transient behaviour starting at the end state of the previous Poincare section. The results of computation the largest Lyapunov exponent presented below were obtained with the use of the classical algorithm based on the differential equations of the tested system. The Gramm-Schmidt reorthonormalization period equal to 0.2 of the period of external excitation was applied. The Lyapunov exponent, which is always zero and related to the external input phase, is omitted here. Figures 2, 3, and 4 contain bifurcation diagrams for variables y 1 and y 2 with increasing bifurcation parameter ξ p in the range 1 . . . 12. In this calculations the nonlocal parameter takes the following values μ = 0, 2.5, 5 nm 2 , respectively. For the case μ = 0 (without consideration small scale effects) one can observe that the coordinate y 1 starting with ξ p = 2.1 passes through the period doubling bifurcations to narrow chaotic zone, while y 2 remains close to zero until ξ p = 6.3, a complication of the y 2 behavior occurs at ξ p > 10.3. Figure 2c shows diagram of the largest Lyapunov exponent associated with the bifurcation diagram presented in Fig. 2a, b. It can be noticed that pos- itive values of the largest Lyapunov exponent confirm domains of chaotic regimes in Fig. 3a, b. Increasing the value of the small-scale parameter μ = 2.5 nm 2 , we observed that the variable y 1 for ξ p > 2.2 and also the variable y 2 for ξ p > 4.1 demonstrate a rich bifurcation scenario, driving to chaotic behaviour by period doubling bifurcation cascade. Furthermore, for ξ p > 8.4 coordinate y 1 tends to zero, whereas y 2 behaves in complicated manner and chaotic zone is alternated by 4periodic window. Mentioned behaviour is confirmed by diagrams of the largest Lyapunov exponent presented in Fig. 3c.
For the higher value of the nonlocal parameter μ = 5 nm 2 one can observe (Fig. 4) a zone of chaotic oscillations, which is preceded by a change in zero solutions and periodic regimes for both coordinates. Moreover, it can be seen the narrow chaotic zone for variable y 1 (ξ p = 7.3), and that y 1 tends to zero for ξ p > 7.6. At the same time, coordinate y 2 undergoes a complex behaviour, where chaotic vibrations are replaced by Fig. 3 Bifurcation diagrams of Poincaré sections (for the coordinates y 1 and y 2 ) and the largest Lyapunov exponent λ 1 with the bifurcation parameter ξ p and for μ = 2.5 nm 2 (isotropic plate) periodic windows and vice versa. Discussed bifurcation scenario is shown in the Fig. 4c with the largest Lyapunov exponent, where zero values correspond to bifurcations of the system.
Next case is aimed at bifurcation analysis with μ as bifurcation parameter in range 0, . . . , 5 nm 2 . Figure 5a, b contains bifurcation diagrams for coordinates y 1 and y 2 with ξ p = 6. One can observe chaotic zone (μ in 0.4..1.5 nm 2 ) for both coordinates with a periodic window. Stating from μ = 3.8 nm 2 variable y 1 becomes close to zero. It can be seen that dynamics reported in these diagrams is confirmed by the largest Lyapunov diagram Fig. 5c, where the chaotic and periodic zones correspond to positive and negative values respectively. For fixed ξ p = 9 rich bifurcation domain can be seen starting from μ = 1.15 nm 2 for coordinate y 1 as well as for coordinate y 2 accompanied by double period bifurcations (Fig. 6a, b). For μ > 1.6 nm 2 variable y 1 tends to zero in contrast to the y 2 variable which has a chaotic behavior with intervals of periodicity. The same Fig. 4 Bifurcation diagrams of Poincaré sections (for the coordinates y 1 and y 2 ) and the largest Lyapunov exponent λ 1 with the bifurcation parameter ξ p and for μ = 5 nm 2 (isotropic plate) chaotic and periodic domains are demonstrated by diagram of the largest Lyapunov exponent, see Fig. 6c. The highest value of the excitation parameter considered is ξ p = 12, Fig. 7a, b. Here we can see two intervals corresponded to chaotic regimes for μ > 0.6 nm 2 and μ > 3.5 nm 2 separated by periodic and narrow chaotic zones. The same scenario is depicted in diagram of the largest Lyapunov exponent, which is presented in Fig. 7c. Thus, we would like to note here that the nonlocal parameter significantly affects the vibration regime.
The appearance of small-scale effects can obviously change the nature of the oscillations with a transition to chaotic behaviour.
The chaotic attractor observed on bifurcation diagrams presented in Fig. 2 without taking into account the small-scale effect for excitation parameter ξ p = 6.2 is considered. The projections of phase plot and Poincaré section are depicted for varable y 1 , see Fig. 8a. Note that y 2 tends to zero for this value of ξ p . The associated largest Lyapunov exponent, as a function of Fig. 5 Bifurcation diagrams of Poincaré sections (for the coordinates y 1 and y 2 ) and the largest Lyapunov exponent λ 1 with the bifurcation parameter μ and for ξ p = 6 (isotropic plate) the number of periods of forcing, is shown in Fig. 8b.  Figures 9a, b and 10a exhibit chaotic attractors corresponding to bifurcation diagrams presented in Fig. 3 (for ξ p = 7), Fig. 4 (for ξ p = 12, the varible y 1 tends to zero for this case), whereas Figs. 9c and 10b demonstrate the corresponding largest Lyapunov exponent. For the chaotic attractor ( μ = 1.2 nm 2 ) detected in bifurcation diagrams with the small scale parameter chosen as bifurcation parameter (see Fig. 5), phase plot, Poincaré section as well as largest Lyapunov exponent are constructed, see Fig. 11. The character of phase plots and Poincaré sections and positive values of the largest Lyapunov exponent presented in Figs. 8,9,10, and 11 represent the chaotic behavior of the small-scale plate for the selected values of the excitation and nonlocal parameters.
Next, we consider the orthotropic graphene sheet with the following mechanical parameters [17,36,37]: For orthotropic plate bifurcation diagrams for coordinates y 1 , y 2 with bifurcation parameter μ (in range 0, . . . , 5 nm 2 ) are constructed taken ξ p = 6 ( Fig. 12), ξ p = 9 (Fig. 13), ξ p = 12 (Fig. 14). Here, as for an isotropic plate, behaviour of an orthotropic graphene sheet is characterized by a complex bifurcation scenario. In this case, a change in the load parameter ξ p significantly affects the location and width of chaotic zones for coordinates y 1 , y 2 . The presented results are confirmed by the corresponded largest Lyapunov exponent diagrams, see Figs. 12c, 13c, 14c. Further, the chaotic attractor (for μ = 1.2 nm 2 ) that appeared on bifurcation diagrams with ξ p = 6 (see Fig. 12) is analysed. The obtained projections of phase plot, Poincaré section for variable y 1 , y 2 as well as the largest Lyapunov exponent, as a function of the number of periods of forcing are presented in Fig. 15. Fig. 7 Bifurcation diagrams of Poincaré sections (for the coordinates y 1 and y 2 ) and the largest Lyapunov exponent λ 1 with the bifurcation parameter μ and for ξ p = 12 (isotropic plate)

Concluding remarks
The presented work is devoted to the study of parametric oscillations of SLGS. The governing equations model geometrically nonlinear vibrations of a smallscale orthotropic plate satisfying the simply supported boundary conditions. The plate subjected to uni-axial in-plane periodic force. The parametric resonance case occurring when the excitation frequency coincides with the doubled natural frequency is studied. The research is based on the use of nonlocal theory, Kirchhoff's hypothesis, and Von Kármán theory. The Bubnov-Galerkin method used in this work in combination with a two-mode approximation of the deflection made it possible to reduce the problem to the analysis of ODEs, the coefficients of which were obtained in an analytical form. By integrating the resulting system of second-order nonlinear equations, the bifurcation analysis was performed. Varying the exciting parameter and the nonlocal parameter, bifurcation diagrams the largest Lyapunov exponent diagrams were constructed, demonstrating a significant qualitative change in the nature of oscillations with an increase in the above mentioned parameters. At the same time, zones of chaotic vibrations were discovered as well as transition scenarios were studied. The use of the size-dependent theory allowed to observe chaotic attractors that were not detected using the not size-dependent classical theory (μ = 0). Phase plots, Poincaré sections and the largest Lyapunov exponent are presented for the purpose of analyzing complex chaotic oscillatory regimes. Poincaré section for the coordinate y 2 (the coordinate y 1 tends to zero) and the largest Lyapunov exponent as a function of number of periods n of parametric forcing for μ = 2.5 nm 2 and ξ p = 12 (isotropic plate) Fig. 11 Phase plots and Poincaré sections for the coordinates y 1 , y 2 and the largest Lyapunov exponent as a function of number of periods n of parametric forcing for μ = 1.2 nm 2 and ξ p = 6 (isotropic plate) Fig. 12 Bifurcation diagrams of Poincaré sections (for the coordinates y 1 and y 2 ) and the largest Lyapunov exponent λ 1 with the bifurcation parameter μ for ξ p = 6 (orthotropic plate) Fig. 13 Bifurcation diagrams of Poincaré sections (for the coordinates y 1 and y 2 ) and the largest Lyapunov exponent λ 1 with the bifurcation parameter μ and for ξ p = 9 (orthotropic plate)

Conflict of interest
The authors declare that they have no conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/ by/4.0/.