Nonlinear vibration of a nonlocal functionally graded beam on fractional visco-Pasternak foundation

This paper investigates the nonlinear dynamic behavior of a nonlocal functionally graded Euler–Bernoulli beam resting on a fractional visco-Pasternak foundation and subjected to harmonic loads. The proposed model captures both, nonlocal parameter considering the elastic stress gradient field and a material length scale parameter considering the strain gradient stress field. Additionally, the von Karman strain–displacement relation is used to describe the nonlinear geometrical beam behavior. The power-law model is utilized to represent the material variations across the thickness direction of the functionally graded beam. The following steps are conducted in this research study. At first, the governing equation of motion is derived using Hamilton’s principle and then reduced to the nonlinear fractional-order differential equation through the single-mode Galerkin approximation. The methodology to determine steady-state amplitude–frequency responses via incremental harmonic balance method and continuation technique is presented. The obtained periodic solutions are verified against the perturbation multiple scales method for the weakly nonlinear case and numerical integration Newmark method in the case of strong nonlinearity. It has been shown that the application of the incremental harmonic balance method in the analysis of nonlocal strain gradient theory-based structures can lead to more reliable studies for strongly nonlinear systems. In the parametric study, it is shown that, on the one hand, parameters of the visco-Pasternak foundation and power-law index remarkable affect the amplitudes responses. On the contrary, the nonlocal and the length-scale parameters are having a small influence on the amplitude–frequency response. Finally, the effects of the fractional derivative order on the system’s damping are displayed at time response diagrams and subsequently discussed.

strain gradient stress field. Additionally, the von Karman strain-displacement relation is used to describe the nonlinear geometrical beam behavior. The powerlaw model is utilized to represent the material variations across the thickness direction of the functionally graded beam. The following steps are conducted in this research study. At first, the governing equation of motion is derived using Hamilton's principle and then reduced to the nonlinear fractional-order differential equation through the single-mode Galerkin approximation. The methodology to determine steady-state amplitude-frequency responses via incremental harmonic balance method and continuation technique is presented. The obtained periodic solutions are verified against the perturbation multiple scales method for the weakly nonlinear case and numerical integration Newmark method in the case of strong nonlinearity. It has been shown that the application of the incremental harmonic balance method in the analysis of nonlocal strain gradient theory-based structures can lead to more reliable studies for strongly nonlinear systems. In the parametric study, it is shown that, on the one hand, parameters of the visco-Pasternak foundation and power-law index remarkable affect the amplitudes responses. On the contrary, the nonlocal and the length-scale parameters are having a small influence on the amplitudefrequency response. Finally, the effects of the fractional derivative order on the system's damping are displayed at time response diagrams and subsequently discussed.

Introduction
Structures with physical properties, which are varied continuously and gradually along a certain direction, are known as functionally graded (FG) materials [18,39,63]. The major advantage of such materials that are emphasized in the engineering practice is that they lack stress concentration, which is a common problem at interfaces of conventional laminated composites [18,39]. FG materials are usually composed of two different material phases such as metal and ceramics. Despite a significant amount of work done in the field of FG structures, there is still a lot of space and need for the investigation of MEMS/NEMS systems based on FG materials. Nanobeams and nanoplates are utilized in different MEMS/NEMS devices [36,49], such as microactuators [21,45], microswitches [73], microsensors [43], nanoscale resonators [20], and energy-harvesting nanodevices [68]. For studying the dynamic behavior of such systems, various approaches based on experiments, molecular dynamics simulations, and continuum mechanics are already employed in the literature [9,15,29,40,42,53]. However, it can be time and skilldemanding to set up and validate the experiment or to implement molecular dynamics simulations for such structures. For these reasons, size-dependent continuum mechanics models gain popularity due to their simplicity in predicting the mechanical behavior of micro-/nano-scaled structural systems. Various nonclassical elasticity theories are used to capture the size effects in micro-and nanostructures. Among them, the most vastly used continuum theories for studying the nanostructures are: nonlocal elasticity theory [48,51,53], strain gradient theory [29], modified couple stress theory (or modified strain gradient theory) [70], and the surface elasticity theory [65]. Some experiments [29] revealed that nonlocal elasticity theory shows limitations in displaying the stiffness-hardening effect. This deficiency can be avoided when nonlocal strain gradient theory (NLSGT) is used as originally described by Lim et al. [37], which includes both nonlocal and length scale effects into the consideration.
Many studies employed NLSGT when analyzing the mechanical behavior of FG structures. Gao et al. [16] investigated the nonlinear free vibration of FG circular nanotubes using NLSGT and two-step perturbation method. Janevski [24,25] studied linear vibration, stability, and buckling of nonlocal strain gradient Euler-Bernoulli and Timoshenko beams under the influence of temperature. El-Borgi et al. [9] investigated the free and forced vibration response of a simply supported FG beam resting on the nonlinear elastic foundation. The authors applied the perturbation method of multiple scales to obtain the amplitude-frequency curves of the system. Other authors [22] studied the heat-induced nonlinear vibration of FG capacitive nanobeam within the framework of NLSGT. The semianalytic perturbation method of averaging was applied to obtain the governing equations and study the steady-state responses. They also used a shooting technique in conjunction with the Floquet theory for capturing the periodic motions and examining their stability. Wang and Shen [66] investigated the lateral nonlinear vibration of an axially moving simply supported viscoelastic nanobeam based on NLSGT. They used a direct multiscale method to obtain the steady-state amplitudefrequency response in the subharmonic parametric resonance state as well as the Routh-Hurwitz criterion to determine the stability of the (non-) zero equilibrium solution. Jalaei et al. [23] investigated the dynamic stability of a temperature-dependent Timoshenko functionally graded nanobeam exposed to the axial excitation load and magnetic field in a thermal environment. The authors used Navier's and Bolotin's method-based approach to solve the problem. Li et al. [33] studied the longitudinal vibration of rods also using the NLSGT and derived analytical solutions for predicting the natural frequencies and mode shapes for specified boundary conditions. They discovered that the NLSGT rod model exerts a stiffness-softening effect when the nonlocal parameter is larger than the length scale parameter and exerts a stiffness-hardening effect in the opposite case. Li [35] investigated the vibration of axially FG beams based on NLSGT and Euler-Bernoulli beam theory and solved the problem via the generalized differential quadrature method. Simsek [59] proposed a beam model for nonlinear free vibration of an FG nanobeam with immovable ends based on the NLSGT and Euler-Bernoulli (EB) beam theory in conjunction with the von Karman's geometric nonlinearity. Liu [38] examined the nonlinear vibrational behavior of FG sandwich NLSGT nanobeams in the presence of initial geometric imperfection. Nonlinearity induced by the von Karman theory and a cosine function similar to the mode shape form is employed to describe the geometric imperfection mode. They used He's variational principle to solve a nonlinear differential equation and obtain nonlinear frequency. Based on NLSGT, Li and Hu [32] and Zhen and Zhou [76] studied the wave propagation in fluid-conveying viscoelastic single-walled carbon nanotubes. Moreover, Li [34] investigated the fluid critical flow velocities of fluid-conveying microtubes modeled using NLSGT and Timoshenko and Euler-Bernoulli beam theories.
One of the pioneering works in the application of fractional calculus in structural mechanics was done by Rossikhin and Shitikova [55]. Their work includes an overview of different papers in this area. The same authors in [54] proposed a methodology based on the Laplace integral transform method to investigate free damped vibrations of diverse linear hereditarily elastic mechanical systems of single and multiple degrees of freedom whose hereditary properties are described by fractional derivatives. Different generalized rheological models were used such as the Maxwell model with one or two fractional parameters (orders of fractional derivatives), the Kelvin-Voigt model, and the standard linear solid model. Later on, Atanackovic and Stankovic studied the existence, regularity, and stability of the solution of an elastic rod on a fractional derivative type of foundation [3] and investigated lateral vibration of the axially loaded rod. Zhang et al. [74] studied the nonlinear dynamic response of a simply supported fractional viscoelastic beam subjected to transverse harmonic excitation. By using the averaging method, the authors obtained a steady-state response of a singlemode system. Numerical results are determined by an algorithm based on the fractional-order Grünwald-Letnikov derivative and verified with analytical results. Eyebe [14] investigated the nonlinear vibration of a nanobeam resting on a fractional-order Winkler-Pasternak foundation by using the D'Alembert principle to obtain the governing equations and a method of multiple scales to approximate the resulting nonlinear problem. Further, Lewandowski [30] investigated the nonlinear, steady-state vibration and stability of harmonically excited fractional viscoelastic beams. The viscoelastic material of the beams is described by using the Zener rheological model with fractional derivatives. Amplitude equations are obtained by using the finite element and the harmonic balance method in conjunction with the continuation method.
In the paper [6], a homogeneous Euler-Bernoulli beam on a Winkler-type nonlinear, viscoelastic and unilateral or bilateral foundation was considered. The presented model was subjected to multiple concentrated or distributed transverse static or dynamic loads. The IHB method was suited for obtaining nonlinear frequency response of the system. Obtained and presented amplitude-frequency diagrams were expectantly similar to diagrams from our study. However, the considered parameters set was different both regarding the loading and the foundation properties. Their interest was the reaction of the foundation, which models foam materials, both to compression and to tension. Instead, we research the influence of the parameters of the fractional-order model of foundation that can represent a range of different materials from foams to rubbers. Furthermore, our numerical calculations were verified with two other numerical methods, namely with multiple time scales and Newmark methods.
In this work, a detailed investigation of the nonlinear vibration of the nonlocal strain gradient Euler-Bernoulli beam resting on the fractional visco-Pasternak foundation and subjected to harmonic loads is performed. The suggested model contains both, nonlocal parameter considering the nonlocal elastic stress field and a material length scale parameter considering the strain gradient stress field. The following steps are conducted in this research study. First, the governing equation of motion is derived using Hamilton's principle and then reduced to the nonlinear fractional-order differential equation via Galerkin approximation. The methodology to determine steady-state amplitude-frequency responses via incremental harmonic balance method and continuation technique is presented. The obtained periodic solutions are verified against the perturbation multiple scales method and numerical integration Newmark method. At last, a detailed parametric study is performed to show the influence of power-law index, nonlocal parameter, length scale parameter, parameters of fractional visco-Pasternak foundation, and load factors on the amplitude-frequency response curves of the proposed nonlinear problem. Additionally, the effects of the fractional derivative order and power-law index on the system's damping are displayed at time response diagrams and subsequently discussed.

Fractional derivative
The vibration of the deformable structures grounded on the different types of foundation is present in a wide range of practical structures. Usually, the impact of the foundation layer has great importance and has to be modeled appropriately. The model of the foundation with different properties can be found in the literature [72]. Visco-Pasternak foundation model used in our study was upgraded with fractional-order time derivatives of the deformation function. This allows us to encompass the whole range of grounds with different properties. We will use the Riemann-Liouville definition (Eq. (1)) when considering the IHB and multiple scales solutions and the Grünwald-Letnikov definition (Eq. (2)) in the case of Newmark method [50,61]. Riemann-Liouville definition is equal to Grünwald-Letnikov definition [50], and these two definitions are equivalent for a wide class of functions and are often used in real physical and engineering problems. For this reason, we can use one definition and then turn to another when calculating the frequency responses by approximate and numerical methods. Here, both definitions are given for clarity.
The left Riemann-Liouville derivative of the continuous and differentiable, on a time interval [a, b], function f (t), is defined as: where α is the fractional-order derivative parameter within the interval 0 < α < 1. Grünwald-Letnikov definition of a fractional derivative is given as: where and [x] means the integer part of x.

Functionally graded material
A small-scale FG beam of width b and thickness h is made of two different materials, and the effective mate- rial properties (e.g., Young's modulus E and density ρ) vary continuously through the beam's thickness (z direction). Those materials' properties regarding geometrical middle axis based on the power-law distribution function of material are: where indices t and b denote the top and bottom layers of the beam and k is the power-law index, which determines the material variation in the thickness direction of the beam. Geometrical and physical middle surfaces of homogeneous materials coincide. However, the change of material properties in one direction shifts the physical middle surface from the geometrical one for some finite length c. Such new system of reference for FG materials and structures is proposed by several authors [4,31,59]. Therefore, to simplify the analysis and avoid the bending-stretching mode coupling effect, we will use a new coordinate system where the x axis lies in the physical middle surface and the vertical axis is given as z, i.e.
Constant c, denoting the position of the physical middle surface, can be calculated as: For our case of a rectangular FG beam with width b and height h, by substituting Eq. (4) in Eq. (6), expression for c is simplified to: By taking the physical middle surface as a reference, material properties can be expressed as:

Nonlocal strain gradient theory
According to the nonlocal strain gradient theory [37], the strain energy U can be expressed as: where σ i j is the nonlocal stress, and σ (1) i jk is the highorder nonlocal stress. Total stress is given by: Constitutive equation for the nonlocal and local parts [12] can be written as: where μ and l are the nonlocal and length scale parameter, respectively, ∇ = ∂ ∂ x , E(z) is the elasticity modulus, ε x x is the axial strain, and ε x x,x is the axial strain gradient. The general nonlocal strain gradient constitutive relation is given as [33]: Salehipour et al. [57] and later Batra [5] proposed the modified nonlocal theory that is applicable to nonhomogenous materials. However, according to [57], when the nabla operator reduces to partial derivative of length coordinate x (∇ = ∂ ∂ x ) and material properties of the beam are only the functions of the thickness coordinate z, then the classical Eringen theory can be used to account for the small-scale effects in FG beams or plates. Moreover, by introducing the physical surface reference system one can avoid coupling between the bending and stretching modes. In our analysis, we adopted both assumptions to study the nonlinear dynamic response of the FG nonlocal beam resting on the fractional visco-Pasternak foundation.

Beam model and equation of motion
Beam model is given in Fig. 2. The displacement field of the Euler-Bernoulli beam is given as: where u x , u y and u z denote the displacements along the length, width and thickness directions, respectively. Terms u and w are the axial and transverse displacements of the physical middle surface, respectively. Thus, the nonzero strain components of Euler-Bernoulli beam with considered geometric nonlinearity take the form: We are considering the following stress resultants: Further, we define the extensional A x x and the bending coefficient D x x as: Note that for homogeneous beam A x x = E A and D x x = E I . By substituting stress resultants Eq. (17) into Eq. (14), the axial force and moment are obtained as: By integrating the general constitutive relation Eq.(14) over area A, or multiplying it with z and integrating over area A, and using relations, Eq.(18) and Eq. (17) lead to The variation of strain energy δU of the FG beam can be given as in [31,60]: Virtual kinetic energy considering both the longitudinal and transverse motions can be given by: In Eq.(23), the mass moments of inertia take the following form: Note that for homogeneous beam m 0 = ρ A and m 2 = ρ I . According to Emam and Nayfeh [10], the first-order mass moment m 1 can be neglected in the virtual kinetic energy (Eq. (23)) since its contribution is relatively small. Virtual work of external loads can be given by [2]: where In Eq. (25), F m is the restoring force due to the visco-Pasternak layer, q is the distributed transverse load, Q is the external shear force, and M is the external bending moment. In Eq.(26), D α is the operator of Riemann-Liouville fractional derivative. In [56], a similar foundation type is introduced but without the fractional time derivatives. Hamilton's principle will be applied by using Eq. (27): By substituting Eqs. (22), (23), and (25) into Eq. (27), the following two equations of motion are obtained: with classical boundary conditions at x = 0 or x = L: and nonclassical boundary conditions at x = 0 or x = L: By assuming the fast dynamics, acceleration in the axial direction in Eq.(28) is negligible. Therefore, N x x = C = const. Substituting Eqs. (14) and (16) into Eq. (17) , the axial force N x x can be written as By substituting Eq.(32) into Eq.(19), one can obtain where and In the case of hinged-hinged beams, the following boundary conditions are valid: Substituting Eq.(34) into Eq.(32) and applying the boundary conditions Eq.(36), one can obtain the expression for the axial force in the following form: Substituting Eq.(37) and the second equation of motion (29) into (21), moment M can be expressed as: Substituting Eq. (38) in (21), we obtain the sizedependent nonlinear equations of motion for an FG Euler-Bernoulli beam model based on the nonlocal strain gradient theory: After substituting relations for the external loads (Eq. (26)) in Eq. (39), it leads to Eq. (40), given by We introduce the following nondimensional parameters: Note that k x , appearing in Eq. (41), is the radius of gyration, defined in Eq. (42) as For the homogenous beam, k x = I x A . Using nondimensional parameters from Eqs. (41) in Eq. (40), nonlinear equation of motion is transformed into the following nondimensional form: The solution of Eq.(43) could be assumed as a sum of products of amplitude and time functions for each mode. The most usual is single-mode discretization, which has been used by many authors (For example [9,14,17,19,22,30,33,37,56,59,64,65]), and solution would be assumed as in Eq. (44). In our case, this is legitimate, since we have only cubic nonlinearity, and Nayfeh and Lacarbonara have shown in their study [46] 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.
where φ n X is the amplitude function, q(τ ) is the time function and n = 1, 2, ... is the mode number. Coefficients s 0 − s 5 are calculated as: By replacing Eq.(44) into Eq.(43) and using Eq. (45), we obtain the following nonlinear fractional-order differential equation where parameters are given as: Note that Eq. (40) can be nondimensionalized in many ways. Among them, the optimal one is given in this paper. Radius of gyration k x (Eq. (42)) is introduced in nondimensionalization process with the purpose to reduce nonlinear parameter θ in Eq. (46). Extreme high values of θ comparing to linear stiffness parameter ω 2 0 would later induce problems with solving fractional-order differential equation of motion (Eq. (46)).

Nonlinear periodic response
Analytical perturbation methods such as the multiple scales method are usually used to solve the nonlinear fractional differential equations in the case of weak nonlinearity [58]. For strong nonlinearities, it is more common to use numerical methods such as the differential quadrature method (DQM) [41] or incremental harmonic balance (IHB) method [47]. A brief review of available numerical methods for solving the aforementioned nonlinear fractional differential equations is given by Zhou et al. [77]. In this study, periodic solutions found by the IHB method are verified with the results from both the perturbation multiple scales and Newmark numerical method. To apply the IHB method, we introduce a new time scale τ = τ into Eq. (46) to obtain the system of nonlinear ordinary differential equations in the following form: For the arbitrarily chosen initial values for q 0 and 0 of the steady-state modal amplitude, a neighboring state of motion has the incremental changes to the current state and can be expressed in the following form: Substituting Eq.(49) into Eq.(48) and neglecting higherorder terms, we obtain a linearized incremental relation given as: where r is residual term given as: To obtain the periodic solutions of the fractional-order differential equation, q 0 and q are expanded as a finite Fourier series of N terms as: [a n cos (nτ ) + b n sin (nτ )] = C A 0 , [ a n cos (nτ ) + b n sin (nτ )] where a 0 a 1 a 2 ... a N b 1 A = [ a 0 a 1 a 2 ... a N b 1 We substitute Eqs. (52), (53), (54), (55) and (56) into Eq. (50) and apply the Galerkin procedure. Since a fractional-order derivative is an aperiodic function, in the integration procedure we choose the time period T = ∞ and average the integration results for the fractional derivative. In the same way, for the periodic function, we choose the time terminal as T = 2π , which leads us to the following system of equations: This gives us a system of linearized algebraic equations in terms of A in the following form: where elements of the Jacobi matrix M, the corrective vector R, and vector V are given in Appendix 1.
In case that we want the solution at a given single frequency, we would set to zero in Eq. (58). Otherwise, we would solve Eq. (58) for both A and , but insert in the first entry of the vector A and transform the system of equations. We initialize solution process by entering guessed values of A and calculate A using Eq. (58). The solution A is then added to the current estimated value of A to determine the new vector A, i.e., A k+1 = A k + A. (59) We repeat this process until the value of the residuum norm |R| is within the preset tolerance (in our case less than 10 −5 ).

The continuation method
For starting the recurrent continuation process, we need to obtain the periodic response in two successive points by using the IHB method. These initial points are usually taken far from the resonant state, where response amplitudes for both of them are having similar and small values. Then, we apply the predictor-corrector method to carry out point-to-point computation for determining the corresponding branches of the frequency responses. Equation (58) can be rewritten in the more general matrix form as: We introduce new vectors X = [ A ] T and X = [ A ] T . Let us also introduce a function g(X ) of vector X in the following form: Eq. (61). Note that the function g(X ) can be defined in many ways, but the one given in Eq. (61) is the most appropriate g(X) = X T X.
(61) We will also introduce arc-length parameter η to follow the direction of the path. An augmented equation would be g(X) − η = 0.
(62) The slope can be determined by using the two previous known points X k−1 and X k−2 on the response curves, such as The first prediction of the next point can be determined by: Equation (60) can be extended with Eq. (62), and then, the tangent stiffness matrix and residual vector can be given in the following form: More information about the continuation method can be found in [6,11,67].

Numerical results
The methodology outlined in the previous section is utilized herein to find the solution of the fractionalorder forced Duffing differential equation Eq. (46) and examine the resonance of a nonlocal FG beam on a fractional visco-Pasternak foundation. The combination of the IHB and path-following methods are introduced to trace branches of periodic solutions of a nonlinear model of a nonlocal strain-gradient beam on a fractional visco-Pasternak foundation with direct transversal harmonic excitation. The obtained diagrams are showing periodic responses given in the form of amplitudefrequency curves. Firstly, beam natural frequency for two different models, obtained by simplifying our model, is verified for the data available in the literature (Tables 1 and 2). In the second part of the numerical study, the validity of the results from the IHB method is examined (Figs. 3, 4, 6), which is then followed by the parametric study in the frequency (Figs. 5,7,8,9,10,11,12,13,14) and time domain (Fig. 15). It is demonstrated that the fractional visco-Pasternak layer has a significant influence on the response amplitudes. Moreover, the results obtained by the IHB method are verified with the results from multiple scales and the Newmark method. The last part of the numerical results section is devoted to the analysis of the influence of different parameters on the response. The results revealed the importance of the first and third harmonics. The parameter values of the presented mechanical model are adopted from the paper [31], extended with parameters for fractional Pasternak layer and FG material, and presented in Table 3. The static part of the excitation force Q 0 is set to zero, and the dynamic part Q 1 is given in the table. When some parameter is varied, remaining coefficients are taken from Table 3. Moreover, it should be noted that the number of adopted harmonics in the Fourier series is N = 6 and this is used in all numerical examples. The amplitudes obtained by the IHB method and corresponding to particular Fourier coefficients Eq. (55) and harmonics Eq. (54) are computed as given in Eq. (66), For verification with the multiple scale method, a small nondimensional bookkeeping parameter takes the value ε = 0.01.    With the purpose of demonstrating the reliability and accuracy of the proposed approach for the determination of the amplitude-frequency responses and corresponding periodic solutions, the obtained results from the IHB are verified with two different approaches-the perturbation method of multiple scales and the direct numerical integration by using the Newmark method. The first one is used to obtain the amplitude-frequency response diagrams, and the second one to capture periodic motions at desired excitation frequencies. The way we applied the Newmark method to solve the nonlinear fractional differential equation of motion Eq. (48) is given in Appendix 3 thoroughly.
First, we will verify the results by comparing the steady-state frequency responses for the superharmonic resonance case = 1 3 ω 0 obtained by the IHB method with the results from the multiple scales method, as given in Figs. 3 and 4. In these figures, response amplitudes corresponding to displacement are given on the ordinate axis, while excitation frequency is on the abscissa.
In Fig. 3, the amplitude-frequency response curves are given for amplitudes A 3 obtained by the IHB, and amplitudes corresponding to the excitation frequency 3 ω 0 obtained by using the multiple scales method. Fractional parameter α is varied. Figure 3b is zoomed Fig. 3a that enables one to clearly compare the obtained results. Data in Fig. 3 reveal that results obtained by these two methods match well. Besides that, we can also observe that an increase of α decreases the amplitude, which is slightly shifted to the right toward higher frequencies.
In Fig. 4, the frequency response curves are given for amplitudes A 3 obtained by using the IHB method and amplitudes corresponding to the excitation frequency = 1 3 ω 0 obtained by using the multiple scales method. External excitation magnitudes are given as: Q 1 = 0.001, Q 1 = 0.002, and Q 1 = 0.003. From this figure, we can observe a good matching between the result obtained by two different methods. Besides that, we can also observe that an increase in the external excitation magnitude increases the amplitude and shift its value to the right toward higher frequencies.
In Fig. 5a and b, the frequency response curves are given for the amplitudes A 1 and A 3 , respectively, which are given on the ordinate axis, while the excitation frequency is on the abscissa. Due to the stiffnesshardening effect of the external excitation force parameter, not only does the maximum amplitude experience a rise but also the frequency response curves are shifted toward higher excitation frequencies. This shifting can be observed for both the first (Fig. 5a)) and third harmonic amplitude (Fig. 5b)). Also, an increase in the external excitation amplitude causes a significant bending of the amplitude-frequency curves so that the multiple-value solutions may exist in the primary resonance case associated with the first and the third harmonic amplitudes. Three periodic orbits are selected from the response curves (marked as star points in Fig. 5), which are then verified with Newmark-based solutions. The periodic solutions are depicted in the phase plane, where the velocity is given on the ordinate axis, while the displacement is given on the abscissa, as shown in sub-figures a, b and c of Fig. 6. We picked two points close and one far from the resonant state (Fig. 5). From Fig. 6, we can observe a good matching between the result obtained by the IHB and Newmark method. However, better overlapping is achieved when we are far from the resonant state.

Parametric study
In the subsequent examples in this chapter, we have shown the influence of different parameters such as nonlocal parameter, strain gradient parameter, powerlaw index, and parameters of fractional visco-Pasternak foundation on amplitude-frequency response. The influence of excitation force is discussed in the previous subsection. Figure 7 shows the amplitude-frequency response of the nonlinear nonlocal strain gradient FG beam on a fractional visco-Pasternak foundation with external excitation for the first A 1 and the third A 3 harmonic amplitudes, and different values of the nonlocal parameter μ. Besides that, selected parts far from and close to resonant state are magnified on figures. From the   Fig. 7, it could be found that variations in the nonlocal parameter are having weak influence in both the first and the third harmonic vibration amplitudes. Due to large nonlinearity and stiffness of the system, influence of the nonlocal parameter on amplitudefrequency response is small. In other words, nonlinearity reduces nonlocal parameter influence on dynamic response of the system.
The amplitude-frequency response curves for different values of the length scale parameter l are given in Fig. 8 for the first A 1 and the third A 3 harmonic amplitudes. Besides that, selected parts far from and close to resonant state are magnified on figures. We observe (a) (b) Fig. 8 The amplitude-frequency response curves of the nonlinear nonlocal strain gradient FG beam on a fractional visco-Pasternak foundation. Amplitudes A 1 (a) and A 3 (b) for different values of the length scale parameter l that variation of the length scale parameter l has a small influence on vibration amplitudes for the primary resonance case and the maximum value. Due to large nonlinearity and stiffness of the system, the influence of the length scale parameter on amplitude-frequency response is small. In other words, nonlinearity reduces length scale parameter influence on dynamic response of the system.
The amplitude-frequency response curves are given for the first A 1 (Fig. 9a)) and the third harmonic amplitude A 3 (Fig. 9b)) and variations of the fractional visco-Pasternak foundation parameter K w . One can observe that an increase of K w decreases the amplitude and (a) (b) Fig. 9 The amplitude-frequency response curves of the nonlinear nonlocal strain gradient FG beam on a fractional visco-Pasternak foundation. Amplitudes A 1 (a) and A 3 (b) are plotted for different values of the parameter K w therefore enlarges the total stiffness of the system. Moreover, an increase of K w as damping parameter decreases the natural frequencies of the system, and therefore, the resonance frequency is shifted to the left. Besides that, by looking at the data in-depth, it can be observed that the angle of curve tilt decreases together with the amplitude toward the curvature center for an increase of K w , which at the same time results in weakening of the hardening-type nonlinear behavior. Figure 10 shows the amplitude-frequency response curves for the first A 1 (Fig. 10a) and the third harmonic amplitude A 3 (Fig. 10b) and variations of the foundation parameter k w . One can observe that an increase of (a) (b) Fig. 10 The amplitude-frequency response curves of the nonlinear nonlocal strain gradient FG beam on a fractional visco-Pasternak foundation. Amplitudes A 1 (a) and A 3 (b) are plotted for different values of the parameter k w k w decreases the amplitude with the stabilizing effect to the system vibrations, and therefore, the total stiffness of the system is enlarged. Besides that, an increase of k w as damping parameter decreases the natural frequencies of the system and shifts the resonance frequency to the right.
The amplitude-frequency response curves in the first A 1 (Fig. 11a) and the third harmonic amplitude A 3 (Fig. 11b) are given for different values of the fractional visco-Pasternak foundation parameter K g . One can notice that an increase of the parameter K g decreases the resonance amplitude that is shifted to the left. This indicates that raise of K g augments the (a) (b) Fig. 11 The amplitude-frequency response curves of the nonlinear nonlocal strain gradient FG beam on a fractional visco-Pasternak foundation. Amplitudes A 1 (a) and A 3 (b) for different values of the parameter K g total stiffness of the system. Furthermore, an increase of the parameter K g causes weakening of the nonlinear hardening behavior of the response. Namely, the hardening-type nonlinearity becomes more apparent when the damping parameter K g is small. Figure 12 shows the amplitude-frequency response curves for the first A 1 (Fig. 12a) and the third harmonic amplitude A 3 (Fig. 12b) for different values of the foundation parameter k g . One can notice that an increase of the parameter k g decreases the resonance amplitude that is shifted to the right significantly enlarging the hardening effects of nonlinearity. This indicates that raise of k g increases the total stiffness of the system.  Fig. 12 The amplitude-frequency response curves of the nonlinear nonlocal strain gradient FG beam on a fractional visco-Pasternak foundation. Amplitudes A 1 (a) and A 3 (b) for different values of the parameter k g By comparing the variation of K w , k w , K g , k g , one can observe that an increase of the parameter K g has a bigger influence on the system's total stiffness increase than the parameter K w , even though both parameters contribute to the amplitude decrease. However, k g and k w have similar effects of moving the amplitudefrequency curve to the right toward the higher values of the external frequency with the light reduction of amplitude values.
The amplitude-frequency response curves in the first A 1 (Fig. 13a) and the third harmonic amplitude A 3 (Fig. 13b)  Pasternak foundation α. It can be noticed that a decrease of the fractional-order parameter α by a step of 0.05 increases the amplitude values by almost double in the primary resonance case. This significant influence of the parameter α is caused by damping features of the system become less pronounced due to the elastic-like behavior of the fractional term. Moreover, a decrease of the fractional derivative parameter α makes the equivalent stiffness coefficient larger, which results in the rightward bending of the amplitude-frequency curves and larger primary resonance frequencies. Amplitudes of the first A 1 (Fig. 14a) and the third harmonic A 3 (Fig. 14b) of the amplitude-frequency response are given for different values of the powerlaw index k that defines the FG material. One can notice that for k = 1 and k = 3 the resonant frequencies and hysteresis domain become larger and more shifted and bent toward the positive direction of lateral axis than for the case when k = 2 and k = 4. This can be attributed to increased stiffness properties of the nonlocal beam for these uneven values of the power-law index that increases the hardening nonlinearity and stiffness features of the system. In this section, we show the time responses of the system obtained via the Newmark method. The influences of the fractional-order derivative parameter α (Fig. 15) are studied to show their effect on the timedependent behavior of the system. To understand the influence of the fractional visco-Pasternak layer on the initial harmonic excitation of the beam, the following initial conditions are adopted q(0) = 1,q(0) = 1. We adopted the following values of fractional parameter α = 0.5, 0.6, 0.7. The dimensionless time period T = 200 is used in this simulation. Similar conclusions can be drawn here as for the amplitude-frequency response. One can observe that an increase of the fractional parameter α leads to stronger damping in time and reduced and attenuated amplitudes of the response. Also, a weak beating phenomenon with decreasing intervals in time can be observed.

Summary of the numerical results
The following conclusions can be drawn about the results presented in this section. The hardening-type nonlinearity becomes more apparent when the force increases and the following parameters decrease: nonlocal parameter μ, strain-gradient length scale parameter l, parameters of visco-Pasternak foundation α, K w , K g . If we increase the external excitation amplitude in this system, the primary resonances will be strength-ened and shifted rightward, i.e., toward higher excitation frequencies. In this case, the hysteresis domain would also increase. Nonlocal and length scale parameters are both having a small influence on the amplitudefrequency response. Parameters of the visco-Pasternak foundation K w and K g augment the total stiffness of the system since their increasing cause response amplitudes decreasing. Specifically, foundation parameters K w and k w have smaller influence on amplitudefrequency response comparing to parameters K g and k g . The even values of power-law index k cause higher amplitude values in comparison with their odd values.
In addition, we remark on the mutual interactions of the regime in the time domain of the single-amplitude mode of beam vibration. This is observed from the amplitude-frequency diagrams of the first mode A 1 ,Figs. 7,8,9,10,11,12,13 and 14, a small jump in amplitude in the region of external frequency about 5, before the resonant region, which corresponds to the contribution of resonant jumps of the third amplitude A 3 . The amplitude of the third mode A 3 has a resonant range around this frequency, and their values go up to values that can be registered for these small jumps on the amplitude diagrams of the first mode A 1 . The second resonant region of the third amplitude A 3 is in the same frequency domain as the one of the first amplitude A 1 , the interval 15-25. Thus, the changes of the A 1 diagram from this interval are contributed also by the behavior of the third amplitude A 3 in this interval.

Conclusions
In this paper, we analyzed the nonlinear vibration problem of a nonlocal beam resting on the fractional visco-Pasternak foundation by using the nonlocal straingradient theory and fractional order damping. The governing equation is derived by using Hamilton's principle and then discretized via the Galerkin approximation, which yields a corresponding nonlinear fractionalorder forced Duffing-type differential equation. The solution is sought for the steady-state superharmonic resonance conditions by using the perturbation multiple time scales method for the weakly nonlinear case and IHB and Newmark method for the strongly nonlinear case. From the verification study, it is revealed that the IHB method is in good agreement with the multiple time scales analysis for the weakly nonlinear case and with the numerical Newmark method for the strongly nonlinear case. The main advantage of the IHB method over the multiple scales method lies in the fact that it does not require an introduction of small parameter, and thus, strong nonlinearity cases can be observed. On the other side, the superiority of the IHB over the Newmark approach is the simple computational implementation and easier determination of periodic solutions. We have also shown that the introduction of the IHB method in the analysis of NLSGT structures can lead to more reliable studies of strongly nonlinear systems. In our parametric study, we concluded that the nonlocal and length scale parameters are having a small influence on the amplitude-frequency response. On the other hand, parameters of the visco-Pasternak foundation remarkably affect the response amplitudes. Finally, the power-law index displays a significant effect on the frequency response, which was also discussed in the numerical analysis. Generally speaking, the system vibration amplitudes are higher for the odd values of the power-law index comparing to materials with the even values of this parameter.
Acknowledgements This research was sponsored by the Serbian Ministry of Education, Science, and Technological Development. D.K. and M.C were funded by the Marie Skłodowska-Curie Actions-European Commission fellowship (Grant No. 799201-METACTIVE, and Grant No. 896942-METASINK, respectively).

Data availability
No datasets are associated with this manuscript. The datasets used for generating the plots and results during the current study can be directly obtained from the numerical simulation of the related mathematical equations in the manuscript.

Conflict of interest:
The authors declare that they have no conflict of interest.

Appendix 1
Elements of the Jacobi matrix M = M 1 + M α 2 , the corrective vector R = R 1 + R α 2 , and vector V = V 1 + V α 2 are defined as: Within each incremental step, only a set of linear equations Eq. (58) has to be solved to obtain the data for the next stage. By applying the procedure established at [47,69]   [ R α 10 = 0, where δ i j is Kronecker delta.

Multiple scales method
Multiple scales is the analytical perturbation method for constructing approximate solutions of nonlinear differential equations. This method is well established in the literature, but it is only valid for small nonlinearities and damping. Therefore, we will use it here only for validation purposes. Equation (46) is well known as the forced Duffing fractional-order differential equation, which can be expressed in terms of small scale parameter as in Eq. (77). Let assume for simplicity Here, we introduce new parameters as γ = γ and θ = θ . The small bookkeeping parameter is put in front of the fractional and nonlinear terms to have weak damping and weak nonlinearity. Please note that the forcing term in Eq. (77) is of the order one (also known as hard forcing), which will help us to study secondary resonances in the system by using the perturbation analysis of the first order. Forcing of order would indicate a primary resonance that is the same as in the Duffing equation [52]. Using the multiple scales method, we will seek the solution of Eq. (77) in the following form: Here, T 0 = τ is the fast time scale and T 1 = τ is the slow time scale. We will analyze the system for superharmonic resonance conditions. Firstly, let us define the time derivatives as d dτ where D n = ∂ ∂ T n , (n = 0, 1, 2, . . . ) and D α−n n+ = ∂ α−n ∂ T α−n n+ , (n = 0, 1, 2, . . .) are classical and Riemann-Liouville's fractional derivative for new time scales [58]. For the fractional derivative of the exponential function [58], restricted to the first-and secondorder approximations, the following relationship will be used: where i is the imaginary unit. Substituting Eqs. (78), (79), (80), (81) into Eq. (77) and then extracting coefficients of 0 and 1 , we obtain the following equations 0 : D 2 0 q 0 + ω 2 0 q 0 = f cos τ, (83) 1 : The solution of Eq. (83) is sought in the form where A is a complex function in terms of slow time scale, and is defined as: Superharmonic resonance 3 ≈ ω 0 Since we have only cubic nonlinearity in Eq. (77), we will consider the case when 3 = ω 0 + σ , where σ is the detuning parameter. By substituting q 0 from Eq. (85) into Eq. (84) and removing the secular terms that grow in time unbounded, i.e., the coefficients of e iω 0 T 0 , we obtain the corresponding solvability conditions as: where A = D 1 A. Then, we use the polar form A = 1 where G L k = (−1) k α k (97) Grunwald-Letnikov coefficients can also be represented in recursive form as: where τ is time step for coarse mesh, and h is time step for fine mesh. Representation of fractional derivative given by Eq.(96) in fine mesh is: where: p is the number of past terms of length h in a time integration step of length τ , j are previous time steps of length τ that can be approximated accurately by a backward Taylor expansion using the displacement, velocity, and acceleration at a certain time step i, k represents overall chunks of j time steps that must be taken into consideration to accurately approximate the fractional derivative at a given point.
Taylor backward expansion for the last j p time steps can be represented as in Eq. (101).
where q i ,q i andq i are displacement, velocity, and acceleration, respectfully, at time step i. Let us neglect higher-order terms. Equation (101) + h −α G L (k−1) j p · · · G L k j p−1 Lets consider equation of motion Eq.(48) in two consecutive time instants. where By substituting Eq.(106) in Eq.(107), we obtain 2 + γ α D 03 q i + γ α D 02 q i where Note that in case of f i = const, Eq.(109) can be solved using the Runge-Kutta method (function ode45 in Matlab). If this is not the case, Eq.(109) can be solved using the Newmark-Beta method. For validation of the IHB solution, the Newmark-Beta method for nonlinear systems is used and implemented according to the procedure presented in [7,13].