11587 FACTA UNIVERSITATIS Series: Mechanical Engineering Vol. 23, No 4, 2025, pp. 709 - 733 https://doi.org/10.22190/FUME230419022N © 2025 by University of Niš, Serbia | Creative Commons License: CC BY-NC-ND Original scientific paper VIBRATION AND STABILITY OF A NONLINEAR NONLOCAL STRAIN-GRADIENT FG BEAM ON A VISCO-PASTERNAK FOUNDATION Nikola Nešić1, Milan Cajić2, Danilo Karličić2, Mihailo Lazarević3, Sondipon Adhikari4 1Faculty of Technical Sciences, University of Priština, Kosovska Mitrovica, Serbia 2Mathematical Institute of the Serbian Academy of Sciences and Arts, Belgrade, Serbia 3Faculty of Mechanical Engineering, University of Belgrade, Serbia 4James Watt School of Engineering, University of Glasgow, United Kingdom ORCID iDs: Nikola Nešić https://orcid.org/0000-0001-6237-4735 Milan Cajić https://orcid.org/0000-0001-5513-0417 Danilo Karličić https://orcid.org/0000-0002-7547-9293 Mihailo Lazarević https://orcid.org/0000-0002-3326-6636 Sondipon Adhikari https://orcid.org/0000-0003-4181-3457 Abstract. This study investigates the stability of periodic solutions of a nonlinear nonlocal strain gradient functionally graded Euler–Bernoulli beam model resting on a visco-Pasternak foundation and subjected to external harmonic excitation. The nonlinearity of the beam arises from the von Kármán strain-displacement relation. Nonlocal stress gradient theory combined with the strain gradient theory is used to describe the stress-strain relation. Variations of material properties across the thickness direction are defined by the power-law model. The governing differential equation of motion is derived by using Hamilton's principle and discretized by the Galerkin approximation. The methodology for obtaining the steady-state amplitude-frequency responses via the incremental harmonic balance method and continuation technique is presented. The obtained periodic solutions are verified against the numerical integration method and stability analysis is performed by utilizing the Floquet theory. Key words: Nonlocal strain gradient theory, Functionally graded beams, Pasternak layer, Duffing oscillator, Incremental harmonic balance method, Floquet multipliers Received: April 19, 2023 / Accepted July 29, 2023 Corresponding author: Nikola Nešić Faculty of Technical Sciences, University of Priština, 38220 Kosovska Mitrovica, Kneza Miloša 7, Serbia E-mail: nikola.nesic@pr.ac.rs 710 N. NEŠIĆ, M. CAJIĆ, D. KARLIČIĆ, M. LAZAREVIĆ, S. ADHIKARI 1. INTRODUCTION Traditional composites are homogeneous mixtures of two or more materials, where compromise is made between the desirable properties of the component materials. On the contrary, a functionally graded material (FGM) is a two-component composite characterized by a compositional gradient from one component to the other and because significant proportions of FGM contain pure forms of each component, the need for compromise is eliminated. In such a case, the properties of both components can be fully utilized [1, 2]. Although a significant amount of work is done in the field of FG structures, there is still a lot of space and need for the investigation of MEMS/NEMS systems composed of FGMs. Nanobeams and nanoplates are used in different MEMS/NEMS devices [3, 4], such as microactuators [5], microswitches [6], micro sensors [7, 8], nanoscale resonators [9], energy harvesting nanodevices [10], etc. Besides that, carbon nanotubes are used as an addition in some mixtures to increase the dynamic stability of the material [11]. The dynamics of such systems can be studied through experiments [12, 13], molecular dynamics simulations [14], and continuum mechanics [15, 16, 17]. The most vastly used continuum theories for studying nanostructures are: nonlocal elasticity theory [17, 18], strain gradient theory [13, 19], modified couple stress theory (or modified strain gradient theory) [20], the surface elasticity theory [21], and the nonlocal strain gradient theory (NLSGT) [22], which includes both nonlocal and length scale effects into consideration. An overview of available theories and methodologies for the analysis of nano-isotropic, nano-functionally graded, and CNT reinforced nanocomposite structures is given by Garg et al. [23]. Nonlinearity can be introduced in beam vibration problems through geometric imperfections in beam material [24] and foundation [25]. Geometric material nonlinearity is usually introduced with von Kármán strain-displacement relation [24, 26, 27]. Simsek [26] examined the nonlinear vibration behavior of an NLSGT Euler-Bernoulli FG beam with von Kármán’s geometric nonlinearity. Liu et al. [24] studied the nonlinear vibration of FG sandwich NLSGT beams in the presence of initial geometric imperfection induced by the von Kármán theory and a cosine function. They solved a nonlinear differential equation and obtained nonlinear frequency relation by using He’s variational principle. A nonlocal theory proposed by Eringen can be utilized in differential and integral forms. The integral form is more general and accurate but complex to apply when solving advanced problems. The differential form is straighter forward in the application. For most cases of boundary conditions, it is scientifically proven that the differential form is equal to the integral form. However, cantilever beams give contradictory results when employing integral and differential forms. Barreta and Marotti de Sciarra [28] suggested a methodology to bridge two forms by replacing the integral form with differential with additional constitutive boundary conditions which is applied to several simple problems [29, 30, 31]. Still, the equation of motion, for a problem similar in complexity to ours, is exclusively derived in differential form as can be found in recently published papers [24, 26, 32, 33, 34]. Therefore, we will utilize differential approach, but limit our research on pinned-pinned, clamped- pinned and clamped-clamped boundary conditions, since clamped-free boundary conditions produce contradictory results. In recent times, stress-driven, strain-driven, and foundation-driven vibrations have become in research focus. Penna et al. [35] investigated nonlinear free vibrations of geometrically imperfect FG nano-beams based on stress-driven nonlocal elasticity. Penna and Feo [36] studied nonlinear free vibrations of functionally graded porous Bernoulli– Vibration and Stability of a Nonlinear Nonlocal Strain-Gradient FG Beam on a Visco-Pasternak... 711 Euler nano-beams resting on a Winkler elastic foundation through a stress-driven nonlocal elasticity model. Vaccaro et al. [37] studied the size-dependent behavior of nonlocal elastic beams by adopting the stress-driven elasticity theory where the kinematics of beams is modeled by the Reddy variational third-order beam theory. Ansari et al. [38] used strain- driven nonlocal formulations of Eringen’s theory in both differential and integral forms for a numerical study of free vibration behavior of piezoelectric Bernoulli–Euler nanoscale beams considering flexoelectric and nonlocal effects. Comparative studies between stress- driven and strain-driven nonlocal elasticity theory have been made by several authors [39, 40, 41]. Elasto-statical displacement-driven nonlocal problems are studied by Vaccaro [42], Patnaik [43] and Pinnola [44]. Moreover, Pinnola et al. [44] have studied the elasto- static problems of classical Euler-Bernoulli beam on a reaction-driven nonlocal foundation. They proposed an integral elasticity mixture considering a convex combination of Winkler (local) and Wieghardt (nonlocal) laws and applied the proposed methodology to find analytical solutions to simple elasto-static problems in structural mechanics. The dynamic behaviour of NLSGT structures composed of FGM was studied by many authors. The most common approach in their analysis includes applying different perturbation methods. For example, Gao et al. [45] studied the nonlinear free vibration of FG circular nanotubes using NLSGT and two steps perturbation method. By using the same methodology, they also studied the vibration of nano arches [46], and FG nanobeams with several different functionally graded distributions [47]. El-Borgi et al. [15] investigated the free and forced vibration response of a simply supported FG beam resting on the nonlinear elastic foundation using the perturbation multiple scales method to obtain the amplitude-frequency curves of the system. Wang and Shen [48] studied the lateral nonlinear vibration of an axially moving simply supported viscoelastic NLSGT beam. For obtaining the steady-state amplitude-frequency responses in the subharmonic parametric resonance regime a direct multiple scales method is used. Jafarsadeghi-Pournaki et al. [49] investigated the heat-induced nonlinear vibration of FG capacitive nanobeam based on NLSGT. They used averaging perturbation method to obtain the governing equations and study the steady-state responses. Li et al. [50] studied analytically the longitudinal vibration of NLSGT rods and discovered that the NLSGT rod model shows a stiffness-softening effect when the nonlocal parameter is larger than the length scale parameter and a stiffness-hardening effect in the opposite case. However, the main disadvantage of perturbation methods is that they can only treat systems with small nonlinearity. Nanostructures modeled by nonlocal stress and strain gradient theory are also studied by other methods, for example, the differential quadrature method. Li [51] used the generalized differential quadrature method to investigate the vibration of axially FG beams based on NLSGT and Euler-Bernoulli beam theory. However, in the recent time, the IHB method becomes popular in the investigation of nonlinear dynamical systems with harmonically varied excitation. In an interesting work of Bhattiprolu et al. [52], they found periodic solutions of a nonlinear Euler-Bernoulli beam resting on viscoelastic uni- and bilateral foundations by using the incremental harmonic balance method (IHBM). However, they investigated vibration of a linear beam on a nonlinear foundation subjected to multiple concentrated forces and only pinned-pinned boundary conditions. Our considered beam is nonlinear with three different boundary conditions. The advantage of IHB over other methods is in its accuracy, the possibility to treat small and huge nonlinearities, and relatively simply implementation. Perturbation methods can solve differential equations with small nonlinear terms and rather simple problems. Therefore, IHB is used in our work to find periodic solutions when strong nonlinearity is present in the system. 712 N. NEŠIĆ, M. CAJIĆ, D. KARLIČIĆ, M. LAZAREVIĆ, S. ADHIKARI Stability of FG nanobeams has been investigated in the literature by using several different criteria and methods. Wang and Shen [48] employed the Routh–Hurwitz criterion to determine the stability of the (non-) zero equilibrium solution of nonlinear laterally vibrating axially moving simply supported viscoelastic NLSGT beam. Jafarsadeghi- Pournaki et al. [49] used a shooting technique in conjunction with the Floquet theory for capturing the periodic motions and examining their stability of heat-induced nonlinear vibration of FG NLSGT capacitive nanobeam. Jalaei et al. [53] investigated the dynamic stability of a temperature-dependent Timoshenko FG nanobeam under the axial excitation load and magnetic field in a thermal environment by using Navier's and Bolotin's method. In this paper, we used Floquet theory to determine the stability of our model because it is suitable to be implemented coupled with the IHB method. Nešić et al. [27] studied the vibration of a nonlocal nonlinear FG beam on the fractional visco-Pasternak foundation based on the amplitude-frequency response analysis and utilizing the IHB, perturbation multiple scales, and Newmark methods to find and confirm the periodic solution. In the current paper, the stability of periodic solutions is investigated for the amplitude-frequency responses of a nonlinear nonlocal strain gradient functionally graded Euler–Bernoulli beam model resting on a visco-Pasternak foundation and subjected to external harmonic excitation. Therefore, this study can be considered as a special case and natural extension of the work done by Nešić et al. in [27] since a special case of the visco-Pasternak foundation model without fractional derivatives is adopted and stability of periodic solutions analysis is performed. Nonlinear von Kármán strain-displacement relation is used together with the nonlocal stress and strain gradient theories. Variations of material properties across the thickness direction are defined by the power-law model. The governing differential equation of motion is derived by using Hamilton's principle and then discretized by the Galerkin approximation. The incremental harmonic balance method and continuation technique are employed to obtain the amplitude-frequency responses and some periodic solutions are verified against the numerical integration Runge-Kutta method. For each periodic solution, a stability check is performed by utilizing Floquet theory. The novelty in presented paper includes introduction of the model of nonlocal nonlinear FG beam vibrating on the visco-Pasternak foundation. Nonlinearity is introduced via von Karman strain-displacement relation and nonlocality through strain gradient and stress gradient constitutive relations. The proposed visco-Pasternak foundation is defined with four parameters enabling modeling of vast range of practical foundations. Originality is obtaining periodic solutions for such a system by using IHB method and continuation technique and analyzing amplitude-frequency response curves for different parameters. Additional novelty is utilizing Floquet theory to determine stability of periodic solutions for different excitation frequencies of such vibrating system and usage of basins of attraction to determine stability of the solution based on initial conditions. Technical interest in considering such structural schemes includes the development of MEMS/NEMS devices composed of small beams that can vibrate on a broad type of elastic supports. In the practical implementation, specific end supports might be preferable. Therefore, three different boundary conditions are considered to demonstrate their influence on vibration behavior. Vibration and Stability of a Nonlinear Nonlocal Strain-Gradient FG Beam on a Visco-Pasternak... 713 2. GOVERNING EQUATION OF FG BEAM ON A VISCO-PASTERNAK FOUNDATION Let us consider a mechanical model of a beam as given in Fig. 1. Relations for Young's modulus and density of FG beam, as well as relations for the nonlocal strain gradient theory can be found in the Appendix of the paper [27]. Fig. 1 Model of the simply-supported nonlinear nonlocal strain gradient FG beam resting on a visco-Pasternak foundation and excited by traverse load The displacement field of the Euler-Bernoulli beam is given as: 𝑢𝑥(𝑥, 𝑧, 𝑡) = 𝑢(𝑥, 𝑡) − 𝑧 𝜕𝑤 𝜕𝑥 , 𝑢𝑦(𝑥, 𝑧, 𝑡) = 0, 𝑢𝑧(𝑥, 𝑧, 𝑡) = 𝑤(𝑥, 𝑡), (1) where ux, uy, and uz denote the displacements along with the length, width, and thickness directions, described by x, y and z coordinates, respectively. Functions u and w are the axial and transverse displacements of the physical middle surface, respectively. The only non- zero strain component of the Euler-Bernoulli beam is described in terms of displacements with von Kármán's nonlinearity as 𝜀𝑥𝑥 = 𝜕𝑢 𝜕𝑥 + 1 2 ( 𝜕𝑤 𝜕𝑥 ) 2 − 𝑧 𝜕2𝑤 𝜕𝑥2 . (2) Following stress resultants are considered in further study: 𝑁𝑥𝑥 = ∫ 𝑡𝑥𝑥𝑑𝐴 𝐴 , 𝑁𝑥𝑥 (0) = ∫ 𝜎𝑥𝑥𝑑𝐴, 𝐴 𝑁𝑥𝑥 (1) = ∫ 𝜎𝑥𝑥 (1) 𝑑𝐴, 𝐴 𝑀 = ∫ 𝑧𝑡𝑥𝑥𝑑𝐴 𝐴 , 𝑀(0) = ∫ 𝑧𝜎𝑥𝑥𝑑𝐴 𝐴 , 𝑀(1) = ∫ 𝑧𝜎𝑥𝑥 (1) 𝑑𝐴. 𝐴 (3) where Nxx, N(0) xx and N(1) xx are axial forces due to total stress txx, nonlocal stress σxx and high-order nonlocal stress σ(1) xx, respectively. All stresses are in axial direction. M, M(0) and M(1) are bending moments resulting from the following stresses in axial direction: total, nonlocal and high-order nonlocal, respectively. Additionally, the extensional Axx and the bending coefficient Dxx are defined in the following form: {𝐴𝑥𝑥 , 𝐷𝑥𝑥} = 𝑏 ∫ {1, 𝑧2} 𝐻/2−𝑐 −𝐻/2−𝑐 𝐸(𝑧)𝑑𝑧, (4) where H represents the beam height and c is the distance between physical and geometrical middle surfaces of the beam. Note that for homogeneous beam Axx = EA and Dxx = EI, with 714 N. NEŠIĆ, M. CAJIĆ, D. KARLIČIĆ, M. LAZAREVIĆ, S. ADHIKARI E denoting Young’s modulus, A cross sectional area and I second moment of area. Reaction force of the visco-Pasternak foundation acting on the beam is defined as [27]: 𝐹𝑚 = (𝑘𝑤 + 𝐾𝑤 𝜕 𝜕𝑡 ) 𝑤 − (𝑘𝑔 + 𝐾𝑔 𝜕 𝜕𝑡 ) 𝜕2𝑤 𝜕𝑥2 , (5) with the elastic Pasternak parameters kw, kg and viscose constants Kw, Kg shown in Fig. 1. Note that the time derivative in Eq. (5) is in the paper [27] replaced with a more general fractional derivative, being a real number, with values between 0 and 1 that have a physical meaning. Foundations described with fractional derivatives as in [27] are more general and therefore can describe a broader range of real value applications. However, a disadvantage of the fractional derivative model is that a stability analysis of fractional differential equations cannot be done with straightforward methods presented in this paper. After following of the procedure presented in [27], we obtain expressions for variations of strain, kinetic and external energy. After applying Hamilton's principle as given in [27], two equations of motion are obtained. We assume fast dynamics, which means that acceleration in an axial direction is negligible. Finally, this leads us to the following equation of motion of a nonlinear FG beam resting on the visco-Pasternak foundation: (1 − 𝜇2 𝜕2 𝜕𝑥2) [− 𝐴𝑥𝑥 2𝐿 ∫ ( 𝜕𝑤 𝜕𝑥 ) 2 𝑑𝑥 𝜕2𝑤 𝜕𝑥2 𝐿 0 − 𝑚2 𝜕4𝑤 𝜕𝑥2𝜕𝑡2 + 𝜕2𝑤 𝜕𝑡2 + 𝑏𝑘𝑤𝑤 + 𝑏𝐾𝑤 𝜕𝑤 𝜕𝑡 −𝑏𝑘𝑔 𝜕2𝑤 𝜕𝑥2 − 𝑏𝐾𝑔 𝜕3𝑤 𝜕𝑡𝜕𝑥2 + 𝑄 𝑐𝑜𝑠 𝛺1 𝑡] + 𝐷𝑥𝑥 (1 − 𝑙2 𝜕2 𝜕𝑥2) 𝜕4𝑤 𝜕𝑥4 = 0, (6) where b is the beam width, μ and l are small parameters in nonlocal and strain gradient theory, m2 is the second moment of inertia, Q and Ω1 are the external force magnitude and frequency, respectively. Eq. (6) is non-dimensional with the following relations: 𝑋 = 𝑥 𝐿 , 𝑊 = 𝑤 𝑘𝑥 , 𝛴 = 𝑙 𝐿 , 𝜆 = 𝜇 𝐿 , 𝜁 = 𝑏 𝐿 , 𝜏 = 𝑡 𝑘𝑥 𝐿2 √ 𝐴𝑥𝑥 𝑚0 , 𝛺 = 𝛺1 𝐿2 𝑘𝑥 √ 𝑚0 𝐴𝑥𝑥 , 𝐹 = 𝑄𝐿4 𝐴𝑥𝑥𝑘𝑥 3 , 𝑦 = 𝑚2 𝑚0𝐿2 , 𝑘𝑥 = √ 𝐷𝑥𝑥 𝐴𝑥𝑥 , 𝑘1 = 𝜁𝑘𝑤𝐿5 𝐴𝑥𝑥𝑘𝑥 2 , 𝐾1 = 𝜁𝐾𝑤𝐿3 𝑘𝑥√𝐴𝑥𝑥𝑚0 , 𝑘2 = 𝜁𝑘𝑔𝐿3 𝐴𝑥𝑥𝑘𝑥 2 , 𝐾2 = 𝜁𝐾𝑔𝐿 𝑘𝑥√𝐴𝑥𝑥𝑚0 . (7) Note that radius of gyration kx, appearing in Eq. (7), for the homogenous beam is kx = (Ix /A)1/2. Non-dimensional form of equation of motion Eq. (6), obtained by substituting Eq. (7) in Eq. (6) is (1 − 𝜆2 𝜕2 𝜕𝑋 2) [− 1 2 ∫ ( 𝜕𝑊 𝜕𝑋 ) 2 𝑑𝑋 𝜕2𝑊 𝜕𝑋 2 1 0 − 𝑦 𝜕4𝑊 𝜕𝑋 2 𝜕𝜏2 + 𝜕2𝑊 𝜕𝜏2 + 𝑘1𝑊 + 𝐾1 𝜕𝑊 𝜕𝜏 −𝑘2 𝜕2𝑊 𝜕𝑋 2 − 𝐾2 𝜕3𝑋 𝜕𝜏𝜕𝑋 2 + 𝐹 𝑐𝑜𝑠 𝛺 𝜏] + (1 − 𝛴2 𝜕2 𝜕𝑋 2) 𝜕4𝑊 𝜕𝑋 4 = 0. (8) The solution of Eq. (8) could be assumed as a sum of products of amplitude and time functions for each mode. Nayfeh and Lacarbonara [54] have shown in their study that in certain cases one-mode Galerkin approximation fails to predict the dynamic behavior of hinged-hinged beams, especially when quadratic type nonlinearity is involved and even modes are observed in certain subharmonic or superharmonic resonance conditions. Since we have only a cubic nonlinear term, we will use only a single mode discretization which Vibration and Stability of a Nonlinear Nonlocal Strain-Gradient FG Beam on a Visco-Pasternak... 715 is the most commonly used method in the investigation of nonlinear structural vibration problems [15, 22, 26, 55, 56]. The solution of Eq. (8) is assumed as 𝑊(𝑋, 𝜏) = 𝜑𝑛(𝑋)𝑞(𝜏), (9) where φn is the amplitude function, q is the time function and n = 1, 2, ... is the mode number. Amplitude function φn, also known as eigenfunction, can be defined as the form of trigonometric functions and should satisfy the boundary conditions. In the present paper, three different boundary conditions are taken into account: pinned-pinned (P-P), clamped- pinned (C-P) and clamped-clamped (C-C), and the amplitude function for each case is given in Eq. (10), respectively [34], as 𝜑𝑛(𝑋) = | | 𝑠𝑖𝑛 𝛾𝑛 𝑋, 𝛾𝑛 = 𝑛𝜋, (𝑠𝑖𝑛 𝛾𝑛 𝑋 − 𝑠𝑖𝑛ℎ 𝛾𝑛 𝑋) − 𝑠𝑖𝑛 𝛾𝑛+𝑠𝑖𝑛ℎ 𝛾𝑛 𝑐𝑜𝑠 𝛾𝑛+𝑐𝑜𝑠ℎ 𝛾𝑛 (𝑐𝑜𝑠 𝛾𝑛 𝑋 − 𝑐𝑜𝑠ℎ 𝛾𝑛 𝑋), 𝑡𝑎𝑛 𝛾𝑛 = 𝑡𝑎𝑛ℎ 𝛾𝑛 , (𝑠𝑖𝑛 𝛾𝑛 𝑋 − 𝑠𝑖𝑛ℎ 𝛾𝑛 𝑋) − 𝑠𝑖𝑛 𝛾𝑛−𝑠𝑖𝑛ℎ 𝛾𝑛 𝑐𝑜𝑠 𝛾𝑛−𝑐𝑜𝑠ℎ 𝛾𝑛 (𝑐𝑜𝑠 𝛾𝑛 𝑋 − 𝑐𝑜𝑠ℎ 𝛾𝑛 𝑋), 𝑐𝑜𝑠 𝛾𝑛 𝑐𝑜𝑠ℎ 𝛾𝑛 = −1. (10) Coefficients s0 - s5 are calculated as {𝑠0, 𝑠1, 𝑠2, 𝑠3, 𝑠4, 𝑠5} = ∫ {𝜑, 𝜑2, 𝜑″𝜑, 𝜑𝐼𝑉𝜑, 𝜑𝑉𝐼𝜑, 𝜑′2}𝑑𝑋. 1 0 (11) We replace Eq. (9) into Eq. (8), and use Eq. (11) to obtain nonlinear differential equation as �̈� + 𝛾�̇� + 𝜔0 2𝑞 + 𝜃𝑞3 = 𝑓 𝑐𝑜𝑠 𝛺 𝜏, (12) with the following coefficients 𝛾 = 𝐾1𝑠1−𝐾2𝑠2−𝜆2𝐾1𝑠2+𝜆2𝐾2𝑠3 𝑠1−𝑦𝑠2−𝜆2𝑠2+𝑦𝜆2𝑠3 , 𝜔0 2 = 𝑘1𝑠1−𝑘2𝑠2−𝜆2𝑘1𝑠2+𝜆2𝑘2𝑠3+𝑠3−𝛴2𝑠4 𝑠1−𝑦𝑠2−𝜆2𝑠2+𝑦𝜆2𝑠3 , 𝜃 = − 1 2 𝑠5𝑠2+ 1 2 𝑠5𝑠3𝜆2 𝑠1−𝑦𝑠2−𝜆2𝑠2+𝑦𝜆2𝑠3 , 𝑓 = −𝑠0𝐹 𝑠1−𝑦𝑠2−𝜆2𝑠2+𝑦𝜆2𝑠3 . (13) 3. NONLINEAR PERIODIC RESPONSE For obtaining the amplitude-frequency response of a nonlinear problem described by Eq. (12), the IHB method and continuation technique will be used. IHB is a vastly used procedure for solving the nonlinear differential equations of Duffing and Mathieu-Duffing type [27, 57]. The continuation technique serves as a supporting methodology for obtaining the periodic solutions in points close to resonant states. 3.1. The incremental harmonic balance method To apply the IHB method, we introduce a new time scale τ ̅ = Ωτ into Eq. (12) to obtain the system of nonlinear ordinary differential equations in the following form 𝛺2 𝑑2𝑞 𝑑𝜏 2 + 𝛾𝛺 𝑑𝑞 𝑑𝜏 + 𝜔0 2𝑞 + 𝜃𝑞3 = 𝑓 𝑐𝑜𝑠 𝜏. (14) 716 N. NEŠIĆ, M. CAJIĆ, D. KARLIČIĆ, M. LAZAREVIĆ, S. ADHIKARI For arbitrarily chosen initial values of q0 and Ω0 for the steady-state modal amplitude, a neighboring state of motion is incrementally changed to the current state and it can be expressed in the following form 𝑞 = 𝑞0 + 𝛥𝑞, 𝛺 = 𝛺0 + 𝛥𝛺. (15) Substituting Eq. (15) into Eq. (14) and neglecting higher-order terms, we obtain a linearized incremental relation given as 𝛺0 2𝛥𝑞″ + 𝛾𝛺0𝛥𝑞′ + 𝜔0 2𝛥𝑞 + 3𝜃𝑞0 2𝛥𝑞 = 𝑟 − 2𝛺0𝑞0 ″ 𝛥𝛺 + 𝑓 𝑐𝑜𝑠 𝜏, (16) where r is the residual term given as 𝑟 = −(𝛺0 2𝑞0 ″ + 𝛾𝛺0𝑞0 ′ + 𝜔0 2𝑞0 + 𝜃𝑞0 3). (17) To obtain the periodic solutions of the differential equation, q0 and Δq are expanded into finite Fourier series of N terms as 𝑞0 = 𝑎0 + ∑ [𝑎𝑛 𝑐𝑜𝑠 𝑛 𝜏 + 𝑏𝑛 𝑠𝑖𝑛 𝑛 𝜏]𝑁 𝑛=1 = 𝑪𝑨𝟎, (18) 𝛥𝑞 = 𝛥𝑎0 + ∑ [𝛥𝑎𝑛 𝑐𝑜𝑠 𝑛 𝜏 + 𝛥𝑏𝑛 𝑠𝑖𝑛 𝑛 𝜏]𝑁 𝑛=1 = 𝑪𝜟𝑨, (19) Where the vector C, and vectors A0 and ΔA with unknown constants are given as: 𝑪 = [1 𝑐𝑜𝑠 𝜏 𝑐𝑜𝑠 2 𝜏 … 𝑐𝑜𝑠 𝑁 𝜏 𝑠𝑖𝑛 𝜏 𝑠𝑖𝑛 2 𝜏 … 𝑠𝑖𝑛 𝑁 𝜏], (20) 𝑨𝟎 = [𝑎0 𝑎1 𝑎2 … 𝑎𝑁 𝑏1 𝑏2 … 𝑏𝑁]𝑇 , (21) 𝜟𝑨 = [𝛥𝑎0 𝛥𝑎1 𝛥𝑎2 … 𝛥𝑎𝑁 𝛥𝑏1 𝛥𝑏2 … 𝛥𝑏𝑁]𝑇 . (22) Substituting Eqs.(18)-(22) into Eq. (16) and applying the Galerkin procedure leads to the system of linearized algebraic equations in terms of ΔA given as 𝑲𝜟𝑨 + 𝑽𝛥𝛺 = 𝑹, (23) where, elements of the Jacobi matrix K of size [2N+1, 2N+1], the corrective vector R, and vector V, both of size [2N+1, 1], are defined as: 𝑲 = 1 2𝜋 ∫ [𝛺0 2𝑪𝑇𝑪″ + 𝛾𝛺0𝑪𝑇𝑪′ + 𝜔0 2𝑪𝑇𝑪 + 3𝜃𝑞0 2𝑪𝑇𝑪] 2𝜋 0 𝑑𝜏, (24) 𝑹 = − 1 2𝜋 ∫ [(𝛺0 2𝑪𝑇𝑪″ + 𝛾𝛺0𝑪𝑇𝑪′ + 𝜔0 2𝑪𝑇𝑪 + 𝜃𝑞0 2𝑪𝑇𝑪)𝑨𝟎 + 𝑪𝑇𝑓 𝑐𝑜𝑠 𝜏] 2𝜋 0 𝑑𝜏, (25) 𝑽 = 1 2𝜋 ∫ [2𝛺𝑪𝑇𝑪″] 2𝜋 0 𝑑𝜏𝑨𝟎, (26) with ( )' denoting the derivative with respect to non-dimensional time τ̅. If we want to get the solution at a given single frequency, we should set ΔΩ to zero in Eq. (23). Otherwise, we solve Eq. (23) for both ΔA and ΔΩ, by inserting ΔΩ in the first entry of the vector ΔA and transforming the system of equations. We initialize the solution process by entering guessed values of A, and calculating ΔA by using the Eq. (23). The solution ΔA is then added to the current estimated value of A to determine the new vector A, i.e, 𝑨𝒌+𝟏 = 𝑨𝒌 + 𝜟𝑨. (27) We repeat this process until the value of the residuum norm |R| is within the preset tolerance (in our case less than 10-10). Vibration and Stability of a Nonlinear Nonlocal Strain-Gradient FG Beam on a Visco-Pasternak... 717 3.1.1. The Continuation Method The main benefit of using the continuation method is the possibility to compute periodic solutions at each point. The requirement for starting the continuation process is the determination of the periodic solutions in two successive points by using the IHB method. It is usual that these initial points are taken far from the resonant state since response amplitudes at both points have similar values, which are also small. Subsequently, the predictor-corrector method can be used to make point-to-point computations and obtain corresponding branches of the amplitude-frequency responses. We introduce new vectors X = [A, Ω]T and ΔX = [ΔA, ΔΩ]T, a function of vector g(X) = XTX, and an arc-length parameter η to follow direction of the path. An augmented equation is given as 𝑔(𝑿) − 𝜂 = 0. (28) When two previous points Xk-1 and Xk-2 of the response curve are known, the slope for computation of the next point can be determined as 𝑿′ = 𝑿𝒌−𝟏−𝑿𝒌−𝟐 ‖𝑿𝒌−𝟏−𝑿𝒌−𝟐‖ . (29) The first prediction of the next point is obtained by 𝑿𝒖 = 𝑿𝒌−𝟏 + 𝛥𝜂𝑿′. (30) After extending Eq. (23) with Eq. (28) we have [ 𝑲 𝑽 𝜕𝑔 𝜕𝑨 𝜕𝑔 𝜕𝛺 ] [ 𝜟𝑨 𝛥𝛺 ] = [ 𝑹 𝛥𝜂 − 𝑔 ]. (31) More information about the continuation method one can find in [52, 58, 59]. 4. STABILITY ANALYSIS When the periodic solution is obtained in the form of Fourier series as in Eq. (18) for a chosen value of Ω, we can investigate the local stability of such periodic solution by using the Floquet theory [60, 61]. Based on this theory, we applied the methodology given by Hsu [62] for the approximation of the transition matrix during one period. Let us consider that the solution of Eq. (14) is given as 𝑞(𝜏) = 𝑞0(𝜏) + 𝛥𝑞 (32) where Δq denotes a small perturbation in the proximity of previously determined periodic solution q0. By substituting Eq. (32) into Eq. (14), and after linearization, one can obtain the linear differential equation with time-periodic coefficients given as: 𝛺2𝛥𝑞″ + 𝛾𝛺𝛥𝑞′ + 𝜔0 2𝛥𝑞 + 3𝜃𝑞0 2𝛥𝑞 = 0, (33) To determine the stability of the obtained periodic solutions using the Floquet theory, the transformation of Eq. (33) into the state-space form should be performed as 𝑑𝒀 𝑑𝜏 = 𝑷(𝜏)𝒀 (34) where Y (τ) = [Δq, Δq']T and P (τ) denotes the periodic matrix with the period T = 2π. 718 N. NEŠIĆ, M. CAJIĆ, D. KARLIČIĆ, M. LAZAREVIĆ, S. ADHIKARI The stability criteria based on the Floquet theory is used for the determination of the local stability of periodic solutions and it is related to the determination of the Floquet multipliers as Hsu gave in [62]. The Floquet multipliers are represented by the eigenvalues of the monodromy matrix. By solving the corresponding eigenvalue problem, for the case when all Floquet multipliers are located inside the unit circle centered at the origin of the complex plane, the determined periodic solutions are stable or asymptotically stable. Otherwise, when the values of Floquet multipliers are outside the unit circle, the periodic solutions are unstable [60]. Depending on where the Floquet multipliers or a pair of complex conjugate multipliers cross the complex plane unit circle, different bifurcation points can be detected. Introduction of the Hsu procedures for the numerical approximation of the transition matrix during one period leads to the numerical determination of Floquet multipliers as a solution of the corresponding eigenvalue problem as Shen et al. stated in [60]. It is assumed that the period T = 2π of the periodic solution q0 (τ) is divided into Nk sub-intervals, where the k-th interval is equal to Δk = τk - τk-1 for τk = kT / N. Furthermore, P (τ) is the continuous periodic matrix with respect to τ and period T, wherein the k-th interval it can be replaced by the constant matrix provided in the case when Nk is chosen to be sufficiently large: 𝑷𝒌 = 1 𝛥𝑘 ∫ 𝑷(𝜏)𝑑𝜏. 𝜏𝑘 𝜏𝑘−1 (35) The monodromy (transition) matrix can be written in the following form 𝑀= ∏ 𝑒𝑷𝒊𝛥𝑖 𝑁𝑘 𝑖=1 ≈ ∏ (𝑰 + ∑ (𝑷𝒊𝛥𝑖)𝑗 𝑗! 𝑁𝑗 𝑗=1 ) 𝑁𝑘 𝑖=1 , (36) where Nj denotes the number of terms in Taylor series, for the approximation of matrix exponent of the constant matrix Pk. From the monodromy matrix M, given in Eq. (36), one can obtain Floquet multipliers as its eigenvalues σ from relation: 𝑑𝑒𝑡(𝑴 − 𝜎𝑰) = 𝟎. (37) For Eq. (14), the periodic matrix P (τ) can be obtained in the following form 𝑷(𝜏) = [ 0 1 − 𝜔0 2+3𝜃𝑞0 2 𝛺2 − 𝛾 𝛺 ]. (38) In this study, for approximation of the monodromy matrix the values Nk = 5000 and Nj = 5 are adopted as the optimal values compromising computational cost and precision. Larger values for Nk and Nj wouldn’t offer some benefit in determination whether point is stable or unstable in this particular numerical example, but would increase computation time significantly. 5. NUMERICAL RESULTS Here, the stability of the obtained periodic solutions of the forced Duffing differential equation Eq. (12) is checked by using the methodology given in the previous section. The IHB and continuation method are numerically applied to trace the periodic solutions branches of a nonlinear model of a nonlocal strain-gradient beam resting on a visco-Pasternak foundation with direct transverse harmonic excitation. Moreover, the Floquet theory is applied to check the stability of periodic solutions for given amplitude-frequency curves. In the first part of the numerical study, we validate our model by comparing the beam's natural frequencies in two Vibration and Stability of a Nonlinear Nonlocal Strain-Gradient FG Beam on a Visco-Pasternak... 719 simplified cases with the results available in the literature (Tables 1 and 2). Then, the validity of the results from the IHB method is examined, which is then followed by the parametric and stability study in the frequency domain. The parameters used in simulations are adopted from the paper [33] and extended with the parameters for the Pasternak layer, FG material, and excitation force as given in Table 3, unless other specified. The number of adopted harmonics in the Fourier series is N = 6. The amplitudes corresponding to particular Fourier coefficients (Eq. (21)) and harmonics (Eq. (20)) are computed as 𝐴0 = 𝑎0, 𝐴𝑖 = √𝑎𝑖 2 + 𝑏𝑖 2, (𝑖 = 1,2, . . . , 𝑁). (39) Table 1 The first five non-dimensional fundamental natural frequencies of a local Euler- Bernoulli beam resting on Winkler-Pasternak foundation for the simply-supported boundary conditions (k1 = 25, k2 = 25) Present Ref.[63] Ref.[64] Ref.[65] Ref.[66] 19.2133 19.2133 19.2133 19.21 19.2178 50.7002 50.7002 50.7002 50.7 50.7804 100.6767 100.677 100.6767 100.7 - 170.0281 170.028 170.0281 170.1 - 258.9868 258.987 258.9868 259.1 Table 2 Comparison of non-dimensional fundamental natural frequencies of simply- supported nonlocal Euler-Bernoulli beam with different nonlocal parameters μ (L = 10, h = 1, ρ = 1, E = 30∙106, ν = 0.3) μ Present Ref.[17] Ref.[67] Ref.[68] 0 9.8293 9.8696 9.8298 9.8696 1 9.3774 9.4159 8.3814 9.4159 2 8.9826 9.0195 8.9892 9.0195 3 8.6338 8.6693 8.6424 8.6693 4 8.3228 8.3569 8.3329 8.3569 Table 3 Parameter values of the presented mechanical model Parameter Symbol Value Unit Young's modulus at top Et 390 GPa Young's modulus at bottom Eb 210 GPa Density at top ρt 3960 kg/m3 Density at bottom ρb 7800 kg/m3 Power-law index k 1 Height of the beam h 100 nm Width of the beam b 1 μm Length of the beam L 10 μm Nonlocal parameter μ 10 nm Length scale parameter l 100 nm Winkler coeff. of viscoelastic layer kw 1e-8 m-1 Winkler coeff. of viscoelastic layer Kw 1e-8 Nsα/m3 Pasternak coeff. of viscoelastic layer kg 1e-8 m Pasternak coeff. of viscoelastic layer Kg 1e-8 Nsα/m Amplitude of excitation force Q 0.002 N 720 N. NEŠIĆ, M. CAJIĆ, D. KARLIČIĆ, M. LAZAREVIĆ, S. ADHIKARI 5.1. Verification Table 1 compares the first five non-dimensional fundamental natural frequencies of a local Euler-Bernoulli beam resting on Winkler-Pasternak foundation for the simply supported boundary conditions with foundation parameters k1 = 25, k2 = 25, taking the values obtained in [63, 64, 65, 66], where fine agreement is achieved. Besides that, we made a comparison of non-dimensional fundamental natural frequencies of simply supported nonlocal Euler-Bernoulli beam with different values of nonlocal stress-gradient parameter μ with data available in the literature [17, 67, 68]. These results are also in good agreement (Table 2). To demonstrate the reliability and accuracy of the presented methodology for the determination of the amplitude-frequency responses and corresponding periodic solutions, the obtained results from the IHB are verified with the direct numerical integration by using the Runge-Kutta method. Choosing the desired IHB periodic solutions at the corresponding excitation frequency (magenta and cyan stars) from Fig. 2 (a, b) and extracting the initial conditions for the use in the Runge-Kutta method, two solutions are compared in Fig. 3 where excellent matching of the results from two different methods IHB and Runge-Kutta can be observed for both excitation frequency examples. Periodic solutions in Fig. 3 are depicted in the phase plane where the velocity is given on the ordinate axis, while the displacement is given on the abscissa. 5.2. Numerical Results of Stability Analysis In this section, we investigate how excitation force, parameters of the functionally graded material, and foundation are affecting the frequency response curves in terms of the first and third amplitudes. For pinned-pinned boundary conditions, Fig. 2 (a, b) shows the frequency response curves for the amplitudes A1 and A3, respectively, which are then given on the ordinate axis while the excitation frequency Ω is on the abscissa. Here, one can observe that the maximum amplitude experiences a rise at two resonance frequencies. Two such regions can be noticed at frequencies around Ω = 4 and Ω = 16-27, depending on Q. In the case of the first dominant amplitude A1, this change is weak in the first resonance region while it is very high for the second region, where the system experiences the hardening stiffness nonlinearity. Here, one can observe jump-down and jump-up behavior as well as coexisting periodic solutions, which is the characteristic of this nonlinear phenomenon [69]. On the other side, despite the much lower values than A1, the amplitudes A3 show an increase in the first resonance region that is more pronounced compared to those occurring within the second region but displaying only the linear-like behavior. However, in the second resonance region, a strong nonlinear hardening stiffness behavior occurs even for very low values of the amplitude. An increase of the excitation force amplitude Q increases the amplitude values and resonant picks. Moreover, a forward frequency sweeping due to an increase of the excitation frequency Ω, starting from small values, shows changes of response amplitudes up to the point when the periodic response loses its stability. This critical point known as fold point is given as upper brown dot in Fig. 2. It can be observed that further forward frequency sweeping produces another branch of stable periodic solutions, which is also known as the lower branch of periodic solutions. Vibration and Stability of a Nonlinear Nonlocal Strain-Gradient FG Beam on a Visco-Pasternak... 721 It can be noticed in Fig. 2 that the upper branch of the amplitude grows until the critical moment (response peak) which represents the maximal value of the response amplitude. Fig. 2 The amplitude-frequency response curves of the NLSGT beam on visco-Pasternak foundation for pinned-pinned (a,b), clamped-pinned (c,d) and clamped-clamped (e,f) boundary conditions. Amplitudes A1 and A3 are plotted for different values of external load parameter Q. The unstable periodic solutions are represented by the tick solid lines. 722 N. NEŠIĆ, M. CAJIĆ, D. KARLIČIĆ, M. LAZAREVIĆ, S. ADHIKARI At this peak point, periodic solutions lose their stability leading to the jump-down phenomena of the response amplitude i.e. a sudden jump of the response amplitude to some lower values. Further forward frequency sweeping produces stable periodic solutions. On the contrary, the backward frequency sweeping process starts from values higher than the resonant frequency. By moving backward, response amplitude increases slowly until the critical point (another fold point) where the periodic solution loses its stability and the jump-up phenomena occur. If we further decrease the excitation frequency another stable (upper) branch of periodic solutions occurs. The branches of unstable periodic solutions displayed as thick solid lines can be detected in between the two different fold points of the same curve by applying the path following method in combination with the Floquet theory. Similar observations can be made for the Fig. 2 (c, f), for clamped-pinned and clamped-clamped boundary conditions. Fig. 3 Periodic response obtained by IHB and RK method for two points Fig. 4 shows the amplitude-frequency responses for four different values of the functional graded material parameter k. For pinned-pinned boundary conditions, the first (Fig. 4 (a)) and the third (Fig. 4 (b)) amplitude are traced, displaying both stable and unstable branches of the response. The linear growth of response amplitude A1 (or A3) appears for an increase of the excitation frequency from some small values until the fold point when the periodic response loses its stability. Same as in the previous case, further forward frequency sweeping produces a lower branch of stable periodic solutions. Therefore, the upper branch of the amplitude grows until the critical moment (response peak) which represents the maximal value of the response amplitude. At this peak point, periodic solutions lose their stability leading to the jump-down phenomena of the response amplitude. Again, the backward frequency sweeping process starting from the values higher than the resonant frequency displays a slow increase until critical point where the periodic solution loses its stability and the jump-up phenomenon occurs. A further decrease of the excitation frequency gives the stable (upper) branch of the periodic responses. The branches of unstable periodic solutions can be detected between two instability points (lower and upper brown points in Fig. 4) and they are represented by the thicker solid lines. In these unstable points, Floquet multipliers exit the unit circle of the complex plane as displayed on the Fig. 7. Similar observations can be made for clamped- pinned and clamped-clamped boundary conditions (Fig. 4 (c, f)). Vibration and Stability of a Nonlinear Nonlocal Strain-Gradient FG Beam on a Visco-Pasternak... 723 Fig. 4 The amplitude-frequency response curves of the nonlinear NLSGT FG beam on visco- Pasternak foundation for pinned-pinned (a,b), clamped-pinned (c,d) and clamped- clamped (e,f) boundary conditions. Amplitudes A1 and A3 are plotted for different values of parameter k. The unstable periodic solutions are represented by the tick solid lines. 724 N. NEŠIĆ, M. CAJIĆ, D. KARLIČIĆ, M. LAZAREVIĆ, S. ADHIKARI Fig. 5 The amplitude-frequency response curves of the nonlinear NLSGT FG beam on visco- Pasternak foundation for pinned-pinned (a,b), clamped-pinned (c,d) and clamped- clamped (e,f) boundary conditions. Amplitudes A1 and A3 are plotted for different values of parameter Kg. The unstable periodic solutions are represented by the tick solid lines. Vibration and Stability of a Nonlinear Nonlocal Strain-Gradient FG Beam on a Visco-Pasternak... 725 Further, Fig. 5 shows amplitude-frequency response curves for changes of the parameter Kg of the visco-Pasternak foundation. One can observe that the change of this parameter yields similar shapes of amplitude-frequency response curves exhibiting the nonlinear hardening stiffness behavior. Moreover, an increase of parameter Kg, related to viscous properties of the foundation, results in a decrease of pick amplitude values. Consequently, this yields shortened unstable branches (thick bold lines) of the periodic solutions located in between the corresponding fold points (Fig. 5). This behavior can be observed for both amplitudes A1 and A3. For clamped-pinned and clamped-clamped boundary conditions, an increase of parameter Kg can even stabilize the system by making an unstable region (bold line) disappear. A more detailed parametric study of amplitude-frequency responses of the similar system with nonlocal strain gradient functionally graded beam on the fractional visco-Pasternak foundation can be found in the paper by Nešić et al. [27]. To study the influence of the type of boundary conditions on the amplitude-frequency response curve three different cases (pinned-pinned, clamped-pinned, and clamped- clamped boundary conditions) are plotted in Fig. 6. It can be observed that reducing degrees of freedom in supports decrease maximal amplitude and shifts resonant frequency to the right, which is expected behavior. Fig. 6 Different boundary conditions In all previous examples, the unstable periodic solutions are detected by the calculation of Floquet multipliers, where at least one of them crosses the unit circle in the complex plane in the +1 direction. Fig. 7 displays a unit circle in the complex plane and Floquet multipliers for the case of parameters adopted from Table 3. It can be observed that Floquet multiplies are crossing the unit circle in the complex plane in the +1 direction for the unstable periodic solutions located between the two different fold points of the amplitude- frequency response curve. It is well known that different bifurcation points can be detected depending on where the Floquet multipliers are crossing the unit circle. Therefore, when Floquet multipliers are crossing the unit circle in the +1 direction, one can expect that a saddle-node bifurcation point occurs in that case. We also observe that number of Floquet multipliers crossing the unit circle is the highest in the case of the simple supports and that the ones that are crossing the unit circle are the most distant to it compared to other 726 N. NEŠIĆ, M. CAJIĆ, D. KARLIČIĆ, M. LAZAREVIĆ, S. ADHIKARI boundary conditions. On the other hand, in the case of clamped-clamped supports, the number of Floquet multipliers crossing the unit circle and their distance to it is minimal. Consequently to this behavior, clamped-clamped boundary conditions contribute to the most and pinned-pinned boundary conditions to the least stable system. Therefore the influence of the boundary conditions on stability cannot be neglected. Fig. 7 Floquet multipliers crossing unit circle (vertical axis is imaginary, horizontal real): a) simply-supported, b) clamped-pinned and c) clamped-clamped boundary conditions Fig. 8 shows a pair of Floquet multipliers for points with different excitation frequency Ω. We can observe in Fig. 8 (a, b) that between approximately Ω = 14.5 and Ω = 23 at least one of two Floquet multipliers has modulus greater than one, meaning that this region is unstable. The stability border is reached when Floquet multipliers modulus (FMM) is equal to one. The degree of instability is defined by the distance of a point in the instability region from the unit circle (Fig. 7) which is equivalent to the distance of slowly growing curve (with small slope) on Fig. 8. For example, for the case of simply-supported beam and Ω = 16 instability is higher than in the case Ω = 22 (Fig. 8 (a, b)). We can also observe from Fig. 8 that with an increase of the system's stiffness by reducing degrees of freedom in boundary conditions, the main instability region is decreased and shifted to the right towards higher excitation frequencies. The stability of the solution for higher excitation frequencies can be achieved by introducing the pinned-pinned boundary conditions. On the contrary, stability on lower excitation frequencies can be reached with clamped-clamped (a) (b) (c) Vibration and Stability of a Nonlinear Nonlocal Strain-Gradient FG Beam on a Visco-Pasternak... 727 boundary conditions. An unstable region with clamped-free boundary conditions is between unstable regions for another two pairs of boundary conditions. This result can be utilized for the practical realization of such a system. Depending on the excitation frequency and operation range of the system, boundary conditions can be selected and physically established so that unstable solutions for amplitudes are avoided, which can prolong the duration of the system. Fig. 8 Floquet multipliers modulus (FMM) of the first (FMM1) and second (FMM2) Floquet multiplier versus excitation frequency Ω Moreover, in the amplitude-frequency responses, one can notice that the coexisting periodic solutions exist for certain frequency ranges. Knowing that the amplitudes of these periodic solutions are dependent on initial conditions and starting positions in excitation frequency sweeping, there are two sets of initial conditions for which the response amplitude will converge towards the lower or upper branch of the periodic solutions forming the basin of attraction as given in Fig. 9 and Fig. 10. Therefore, the system will possess two attractors for the given excitation frequency and system parameters. These two attractors are highlighted in two different colors for two sets of initial conditions. For the presented nonlinear system, the basin of attraction is given in Fig. 9 and Fig. 10 and for the excitation frequency Ω = 20 and Ω = 25 and coexisting periodic solutions from the 728 N. NEŠIĆ, M. CAJIĆ, D. KARLIČIĆ, M. LAZAREVIĆ, S. ADHIKARI amplitude-frequency response curve given in Fig. 2 and varying initial conditions. One can observe that the basin of attraction is displaying two colors with two coexisting periodic solutions in Fig. 9 and Fig. 10. Here, the case with the initial condition highlighted in red color converges to the upper stable branch while the case with the initial conditions highlighted in blue converges to the lower stable branch. Since excitation frequency is close to the bifurcation point, the solution converges to both upper and lower stable branches of periodic solutions. In the special case of simple-supported beam, the set of initial conditions that converges to the upper stable branch of periodic solutions is dominating over the set of initial conditions that converges to the lower branch of solutions. Fig. 9 Basin of attraction for Ω = 20 Fig. 10 Basin of attraction for Ω = 25 for clamped-clamped boundary conditions To summarize, in this chapter, after a short verification of eigenfrequencies for simplified model with the results available in the literature, stability analysis is performed. For different values of excitation frequency acting in transversal direction, the first (A1) and the third (A3) harmonic of the amplitudes for displacement of the beam in transversal direction are obtained and given in the Figs. 2, 4, 5. These figures, for given frequency, give us information about the number of solutions (1, 2 or 3) for each of the amplitudes Vibration and Stability of a Nonlinear Nonlocal Strain-Gradient FG Beam on a Visco-Pasternak... 729 (A1 and A3) and whether the solution is stable or not. All this is to predict the vibration behavior of the beam at certain frequency and set of the parameters. Selected solution for amplitude is stable if a pair of Floquet multipliers is in the unit circle of the complex plane as presented on the Fig. 7. Floquet theory is basically a tool for proving stability. In Figs 7-8 are graphically given values of Floquet multipliers for different values of excitation frequency from the range Ω = 0 – 30, enabling reader to determine how far from stability border are solutions at particular excitation frequency. Finally, for specific frequency, which is in the range of multiple amplitude solutions, basins of attractions are given in Figs. 9-10, describing whether amplitude will be greater or smaller, depending on initial conditions. 6. CONCLUSIONS In this work, a stability analysis of a nonlinear nonlocal strain-gradient beam resting on the visco-Pasternak foundation is performed. The governing equation of the system is derived by using Hamilton's principle and then discretized by using the Galerkin approximation method. The periodic solutions of the strongly nonlinear differential equation are then sought by using the incremental harmonic balance method in combination with the continuation technique. The results obtained by the proposed approach are verified with the direct numerical integration technique (Runge-Kutta method) and good agreement is confirmed. A parametric study is performed through the amplitude-frequency responses and the stability of each periodic solution is examined by using the Floquet stability theory and Hsu procedures. The main contributions of this work can be summarized as follows: ▪ The parametric study is performed and stability of obtained periodic solutions of a nonlocal strain-gradient functionally graded beam on the visco-Pasternak foundation is investigated. It is demonstrated that the hardening stiffness nonlinearity becomes more apparent for an increase of the excitation amplitude and for lower values of the parameter of functionally graded material, which consequently leads to an increase of the unstable branches of periodic solutions. Opposite to this, an increase of the parameter related to viscous properties of the foundation yields lower picks in the amplitude-frequency response and smaller unstable branches, therefore increasing the stability of the system. ▪ Investigation of the Floquet multipliers and their behavior in the complex plane revealed the nature of bifurcation points where unstable solution appears. Moreover, the exiting of Floquet multipliers out of the unit circle in the +1 direction demonstrated the existence of saddle-node bifurcation points. ▪ Modulus of Floquet multipliers is compared for different values of excitation frequencies and shown that in an unstable region, instability degree is higher for lower excitation frequencies. ▪ To obtain stable solutions on higher excitation frequencies, pinned-pinned boundary conditions should be applied. On the contrary, stability on lower excitation frequencies can be reached with clamped-clamped boundary conditions. An unstable region with clamped-free boundary conditions is between unstable regions for another two pairs of boundary conditions. ▪ The occurrence of coexisting periodic solutions and hardening stiffness effect in the strongly nonlinear functionally graded beam system is confirmed as well as their dependence on the initial conditions and starting positions of the excitation frequency sweeping. Therefore, the basin of attraction is obtained for such a system to examine which of the stable periodic solutions will prevail for different sets of initial conditions. 730 N. NEŠIĆ, M. CAJIĆ, D. KARLIČIĆ, M. LAZAREVIĆ, S. ADHIKARI The paper illustrates how the incremental harmonic balance method in combination with the continuation technique and Floquet theory can be efficiently utilized to study the stability of strongly nonlinear structural vibration systems with different sizes, boundary conditions, and foundation effects included. The future perspectives of the proposed research and the developed model is application of the model and methodology presented in this paper in the development of MEMS devices, since the presented model describes oscillating nonlocal beam, and nonlocal theory is used in small scale systems. Additionally, when neglecting two nonlocal parameters, presented model of beam oscillating on visco- Pasternak foundation with four adjustable parameters can be used to describe wide range of foundations used in civil and mechanical engineering and oscillatory behavior of the nonlinear beam resting on it. Acknowledgement: Authors N. Nešić, D. Karličić and M. Lazarević acknowledge the support by the Serbian Ministry of Education, Science, and Technological Development through the Faculty of Technical Sciences, University of Priština in Kosovska Mitrovica, the Mathematical Institute of the Serbian Academy of Sciences and Art and Faculty of Mechanical Engineering, Univerisity of Belgrade (grant no. 451-03-68/2022-14/200105). Authors M. Cajić and S. Adhikari acknowledge funding from European Union’s Horizon 2020 research and innovation program under the Marie Skłodowska-Curie grant agreement No. 896942 (METASINK). REFERENCES 1. Ghayesh, M.H., Farajpour, A., 2019, A review on the mechanics of functionally graded nanoscale and microscale structures, International Journal of Engineering Science, 137, pp. 8-36. 2. Mahamood, R.M., Akinlabi, E.T., Shukla, M., Pityana, S.L., 2012, Functionally graded material: an overview, Proceedings of the World Congress on Engineering, Vol. 3. London, UK: International Association of Engineers (IAENG). 3. Peddieson, J., Buchanan, G.R., McNitt, R.P., 2003, Application of nonlocal continuum models to nanotechnology, International journal of engineering science, 41(3-5), pp. 305-312. 4. Li, X., Bhushan, B., Takashima, K., Baek, C.W., Kim, Y.K., 2003, Mechanical characterization of micro/nanoscale structures for MEMS/NEMS applications using nanoindentation techniques, Ultramicroscopy, 97(1-4), pp. 481-494. 5. Najar, F., Nayfeh, A.H., Abdel-Rahman, E.M., Choura, S., El-Borgi, S., 2010, Dynamics and global stability of beam-based electrostatic microactuators, Journal of Vibration and Control, 16(5), pp. 721-748. 6. Zavracky, P.M., McGruer, N.E., Morrison, R.H., Potter, D., 1999, Microswitches and microrelays with a view toward microwave applications, International Journal of RF and Microwave Computer‐Aided Engineering, 9(4), pp. 338-347. 7. Moser, Y., Gijs, M.A., 2007, Miniaturized flexible temperature sensor, Journal of Microelectromechanical Systems, 16(6), pp. 1349-1354. 8. Ceballes, S., Abdelkefi, A., 2021, Uncertainty analysis and stochastic characterization of carbon nanotube-based mass sensor with multiple deposited nanoparticles, Sensors and Actuators A: Physical, 332(2), 113182. 9. He, X.Q., Kitipornchai, S., Liew, K.M., 2005, Resonance analysis of multi-layered graphene sheets used as nanoscale resonators, Nanotechnology, 16(10), 2086. 10. Wang, Z.L., Wu, W., 2012, Nanotechnology‐enabled energy harvesting for self‐powered micro‐ /nanosystems, Angewandte Chemie International Edition, 51(47), pp. 11700-11721. 11. Khan, A.A., Khan, I.U., Noman, M., Khan, U.A. and Bilal, M., 2022, Performance Evaluation of Flexible Pavement Using Carbon Nanotubes and Plastic Waste as Admixtures, Tehnički vjesnik, 29(1), pp. 9-14. 12. Moory-Shirbani, M., Sedighi, H.M., Ouakad, H.M., Najar, F., 2018, Experimental and mathematical analysis of a piezoelectrically actuated multilayered imperfect microbeam subjected to applied electric potential, Composite Structures, 184, pp. 950-960. 13. Lam, D.C., Yang, F., Chong, A.C.M., Wang, J., Tong, P., 2003, Experiments and theory in strain gradient elasticity, Journal of the Mechanics and Physics of Solids, 51(8), pp. 1477-1508. Vibration and Stability of a Nonlinear Nonlocal Strain-Gradient FG Beam on a Visco-Pasternak... 731 14. Mehralian, F., Beni, Y.T., Zeverdejani, M.K., 2017, Calibration of nonlocal strain gradient shell model for buckling analysis of nanotubes using molecular dynamics simulations, Physica B: Condensed Matter, 521, pp. 102-111. 15. El-Borgi, S., Fernandes, R., Reddy, J.N., 2015, Non-local free and forced vibrations of graded nanobeams resting on a non-linear elastic foundation, International Journal of Non-Linear Mechanics, 77, pp. 348-363. 16. Fan, Y., Xiang, Y., Shen, H.S., 2020, Nonlinear dynamics of temperature-dependent FG-GRC laminated beams resting on visco-Pasternak foundations, International Journal of Structural Stability and Dynamics, 20(01), 2050012. 17. Reddy, J., 2007, Nonlocal theories for bending, buckling and vibration of beams, International journal of engineering science, 45(2-8), pp. 288-307. 18. Pradhan, S.C., Phadikar, J.K., 2009, Nonlocal elasticity theory for vibration of nanoplates, Journal of Sound and Vibration, 325(1-2), pp. 206-223. 19. Limkatanyu, S., Sae-Long, W., Rungamornrat, J., Buachart, C., Sukontasukkul, P., Keawsawasvong, S. and Chindaprasirt, P., 2022. Bending, buckling and free vibration analyses of nanobeam-substrate medium systems, Facta Universitatis, Series: Mechanical Engineering, 20(3), pp. 561-587. 20. Yang, F.A.C.M., Chong, A.C.M., Lam, D.C.C., Tong, P., 2002, Couple stress based strain gradient theory for elasticity, International journal of solids and structures, 39(10), pp. 2731-2743. 21. Wang, G.F., Feng, X.Q., 2007, Effects of surface elasticity and residual surface tension on the natural frequency of microbeams, Applied physics letters, 90(23), 231904. 22. Lim, C.W., Zhang, G., Reddy, J., 2015, A higher-order nonlocal elasticity and strain gradient theory and its applications in wave propagation, Journal of the Mechanics and Physics of Solids, 78, pp. 298-313. 23. Garg, A., Chalak, H.D., Zenkour, A.M., Belarbi, M.O., Houari, M.S.A., 2022, A review of available theories and methodologies for the analysis of nano isotropic, nano functionally graded, and CNT reinforced nanocomposite structures, Archives of Computational Methods in Engineering, 29, pp. 2237–2270. 24. Liu, H., Lv, Z., Wu, H., 2019, Nonlinear free vibration of geometrically imperfect functionally graded sandwich nanobeams based on nonlocal strain gradient theory, Composite Structures, 214, pp. 47-61. 25. Fallah, A., Aghdam, M.M., 2011, Nonlinear free vibration and post-buckling analysis of functionally graded beams on nonlinear elastic foundation, European Journal of Mechanics-A/Solids, 30(4), pp. 571-583. 26. Şimşek, M., 2016, Nonlinear free vibration of a functionally graded nanobeam using nonlocal strain gradient theory and a novel Hamiltonian approach, International Journal of Engineering Science, 105, pp. 12-27. 27. Nešić, N., Cajić, M., Karličić, D., Obradović, A., Simonović, J., 2022, Nonlinear vibration of a nonlocal functionally graded beam on fractional visco-Pasternak foundation, Nonlinear Dynamics, 107(3), pp. 2003-2026. 28. Barretta, R., de Sciarra, F.M., 2018, Constitutive boundary conditions for nonlocal strain gradient elastic nano- beams, International Journal of Engineering Science, 130, pp. 187-198. 29. Apuzzo, A., Barretta, R., Faghidian, S.A., Luciano, R., De Sciarra, F.M., 2019, Nonlocal strain gradient exact solutions for functionally graded inflected nano-beams, Composites Part B: Engineering, 164, pp. 667-674. 30. Barretta, R., Faghidian, S.A., Luciano, R., Medaglia, C.M., Penna, R., 2018, Free vibrations of FG elastic Timoshenko nano-beams by strain gradient and stress-driven nonlocal models, Composites Part B: Engineering, 154, pp. 20-32. 31. Pinnola, F.P., Faghidian, S.A., Barretta, R., de Sciarra, F.M., 2020, Variationally consistent dynamics of nonlocal gradient elastic beams, International Journal of Engineering Science, 149, 103220. 32. Sahmani, S., Safaei, B., 2020, Influence of homogenization models on size-dependent nonlinear bending and postbuckling of bi-directional functionally graded micro/nano-beams, Applied Mathematical Modelling, 82, pp. 336-358. 33. Li, L., Hu, Y., 2016, Nonlinear bending and free vibration analyses of nonlocal strain gradient beams made of functionally graded material, International Journal of Engineering Science, 107, pp. 77-97. 34. Liu, H., Wu, H., Lyu, Z., 2020, Nonlinear resonance of FG multilayer beam-type nanocomposites: effects of graphene nanoplatelet-reinforcement and geometric imperfection, Aerospace Science and Technology, 98, 105702. 35. Penna, R., Feo, L., Fortunato, A., Luciano, R., 2021, Nonlinear free vibrations analysis of geometrically imperfect FG nano-beams based on stress-driven nonlocal elasticity with initial pretension force, Composite Structures, 255, 112856. 36. Penna, R., Feo, L., 2020, Nonlinear dynamic behavior of porous and imperfect Bernoulli-Euler functionally graded nanobeams resting on Winkler elastic foundation, Technologies, 8(4), 56. 37. Vaccaro, M.S., Barretta, R., Marotti de Sciarra, F., Reddy, J.N., 2022, Nonlocal integral elasticity for third-order small-scale beams, Acta Mechanica, 233(6). pp. 2393–2403. 38. Ansari, R., Faraji Oskouie, M., Nesarhosseini, S., Rouhi, H., 2022, Vibrations of piezoelectric nanobeams considering flexoelectricity influence: a numerical approach based on strain-driven nonlocal differential/integral models, Journal of the Brazilian Society of Mechanical Sciences and Engineering, 44(2), 57. 732 N. NEŠIĆ, M. CAJIĆ, D. KARLIČIĆ, M. LAZAREVIĆ, S. ADHIKARI 39. Romano, G., Barretta, R., 2017, Stress-driven versus strain-driven nonlocal integral model for elastic nano- beams, Composites Part B: Engineering, 114, pp. 184-188. 40. Faraji Oskouie, M., Ansari, R., Rouhi, H., 2018, Bending of Euler–Bernoulli nanobeams based on the strain-driven and stress-driven nonlocal integral models: a numerical approach, Acta Mechanica Sinica, 34(5), pp. 871-882. 41. Faraji Oskouie, M., Ansari, R., Rouhi, H., 2018, A numerical study on the buckling and vibration of nanobeams based on the strain and stress-driven nonlocal integral models, International Journal of Computational Materials Science and Engineering, 7(03), 1850016. 42. Vaccaro, M.S., Pinnola, F.P., Marotti de Sciarra, F., Barretta, R., 2021, Elastostatics of Bernoulli–Euler beams resting on displacement-driven nonlocal foundation, Nanomaterials, 11(3), 573. 43. Patnaik, S., Sidhardh, S., Semperlotti, F., 2022, Displacement-driven approach to nonlocal elasticity, European Journal of Mechanics-A/Solids, 92, 104434. 44. Pinnola, F.P., Vaccaro, M.S., Barretta, R., Marotti de Sciarra, F., Ruta, G., 2023, Elasticity problems of beams on reaction-driven nonlocal foundation, Archive of Applied Mechanics, 93, pp. 41–71. 45. Gao, Y., Xiao, W., Zhu, H., 2019, Nonlinear vibration of functionally graded nano-tubes using nonlocal strain gradient theory and a two-steps perturbation method, Structural Engineering and Mechanics, 69(2), pp. 205-219. 46. Gao, Y., Xiao, W.S., Zhu, H., 2020, Snap-buckling of functionally graded multilayer graphene platelet- reinforced composite curved nanobeams with geometrical imperfections, European Journal of Mechanics- A/Solids, 82, 103993. 47. Gao, Y., Xiao, W.S., Zhu, H., 2019, Nonlinear vibration analysis of different types of functionally graded beams using nonlocal strain gradient theory and a two-step perturbation method, The European Physical Journal Plus, 134(1), 23. 48. Wang, J., Shen, H., 2019, Nonlinear vibrations of axially moving simply supported viscoelastic nanobeams based on nonlocal strain gradient theory. Journal of Physics: Condensed Matter, 31(48), 485403. 49. Jafarsadeghi-Pournaki, I., Azizi, S., Zamanzadeh, M., Madinei, H., Shabani, R., Rezazadeh, G., 2020, Size- dependent dynamics of a FG nanobeam near nonlinear resonances induced by heat, Applied Mathematical Modelling, 86, pp. 349-367. 50. Li, L., Hu, Y., Li, X., 2016, Longitudinal vibration of size-dependent rods via nonlocal strain gradient theory, International Journal of Mechanical Sciences, 115-116, pp. 135-144. 51. Li, X., Li, L., Hu, Y., Ding, Z., Deng, W., 2017. Bending, buckling and vibration of axially functionally graded beams based on nonlocal strain gradient theory, Composite Structures, 165, pp. 250-265. 52. Bhattiprolu, U., Bajaj, A.K., Davies, P., 2016, Periodic response predictions of beams on nonlinear and viscoelastic unilateral foundations using incremental harmonic balance method, International Journal of Solids and Structures, 99, pp. 28-39. 53. Jalaei, M.H., Arani, A.G., Nguyen-Xuan, H., 2019, Investigation of thermal and magnetic field effects on the dynamic instability of FG Timoshenko nanobeam employing nonlocal strain gradient theory, International Journal of Mechanical Sciences, 161-162, 105043. 54. Nayfeh, A.H., Lacarbonara, W., 1997, On the discretization of distributed-parameter systems with quadratic and cubic nonlinearities, Nonlinear Dynamics, 13(3), pp. 203-220. 55. Hashemian, M., Falsafioon, M., Pirmoradian, M., Toghraie, D., 2020, Nonlocal dynamic stability analysis of a Timoshenko nanobeam subjected to a sequence of moving nanoparticles considering surface effects, Mechanics of Materials, 148, 103452. 56. Trabelssi, M., El-Borgi, S., Friswell, M.I., 2020, A high-order FEM formulation for free and forced vibration analysis of a nonlocal nonlinear graded Timoshenko nanobeam based on the weak form quadrature element method, Archive of Applied Mechanics, 90(10), pp. 2133-2156. 57. Nešić, N., Cajić, M., Karličić, D., Janevski, G., 2021, Nonlinear superharmonic resonance analysis of a nonlocal beam on a fractional visco-Pasternak foundation, Proceedings of the Institution of Mechanical Engineers, Part C: Journal of Mechanical Engineering Science, 235(20), pp. 4594-4611. 58. Engelnkemper, S., Gurevich, S.V., Uecker, H., Wetzel, D., Thiele, U., 2019, Continuation for thin film hydrodynamics and related scalar problems, Computational modelling of bifurcations and instabilities in fluid dynamics, Computational Methods in Applied Sciences, 50, pp. 459-501 59. Wang, S., Hua, L., Yang, C., Han, X.. Su, Z., 2019, Applications of incremental harmonic balance method combined with equivalent piecewise linearization on vibrations of nonlinear stiffness systems, Journal of Sound and Vibration, 441, pp. 111-125. 60. Shen, J.H., Lin, K.C., Chen, S.H., Sze, K.Y., 2008, Bifurcation and route-to-chaos analyses for Mathieu–Duffing oscillator by the incremental harmonic balance method, Nonlinear Dynamics, 52(4), pp. 403-414. 61. Huang, J.L., Su, R.K.L., Lee, Y.Y., Chen, S.H., 2011, Nonlinear vibration of a curved beam under uniform base harmonic excitation with quadratic and cubic nonlinearities, Journal of Sound and Vibration, 330(21), pp. 5151-5164. Vibration and Stability of a Nonlinear Nonlocal Strain-Gradient FG Beam on a Visco-Pasternak... 733 62. Hsu, C.S., 1974, On approximating a general linear periodic system, Journal of Mathematical Analysis and Applications, 45(1), pp. 234-251. 63. Togun, N., Bağdatlı, S.M., 2016, Nonlinear vibration of a nanobeam on a Pasternak elastic foundation based on non-local Euler-Bernoulli beam theory, Mathematical and Computational Applications, 21(1), 3. 64. Eyebe, G.J., Betchewe, G., Mohamadou, A., Kofane, T.C., 2018, Nonlinear vibration of a nonlocal nanobeam resting on fractional-order viscoelastic Pasternak foundations, Fractal and Fractional, 2(3), 21. 65. Yokoyama, T., 1987, Vibrations and transient responses of Timoshenko beams resting on elastic foundations, Ingenieur-Archiv, 57(2), pp. 81-90. 66. Mustapha, K.B., Zhong, Z.W., 2010, Free transverse vibration of an axially loaded non-prismatic single-walled carbon nanotube embedded in a two-parameter elastic medium, Computational Materials Science, 50(2), pp. 742-751. 67. Ansari, R., Gholami, R., Hosseini, K., Sahmani, S., 2011, A sixth-order compact finite difference method for vibrational analysis of nanobeams embedded in an elastic medium based on nonlocal beam theory, Mathematical and Computer Modelling, 54(11-12), pp. 2577-2586. 68. Zhao, X., Wang, C.F., Zhu, W.D., Li, Y.H., Wan, X.S., 2021, Coupled thermoelastic nonlocal forced vibration of an axially moving micro/nano-beam, International Journal of Mechanical Sciences, 206, 106600. 69. Brennan, M.J., Kovacic, I., Carrella, A., Waters, T.P., 2008, On the jump-up and jump-down frequencies of the Duffing oscillator, Journal of Sound and Vibration, 318(4-5), pp. 1250-1261.