Dynamic stability of a nonlinear multiple-nanobeam system

We use the incremental harmonic balance (IHB) method to analyse the dynamic stability problem of a nonlinear multiple-nanobeam system (MNBS) within the framework of Eringen’s nonlocal elasticity theory. The nonlinear dynamic system under consideration includes MNBS embedded in a viscoelastic medium as clamped chain system, where every nanobeam in the system is subjected to time-dependent axial loads. By assuming the von Karman type of geometric nonlinearity, a system of m nonlinear partial differential equations of motion is derived based on the Euler–Bernoulli beam theory and D’ Alembert’s principle. All nanobeams in MNBS are considered with simply supported boundary conditions. Semi-analytical solutions for time response functions of the nonlinear MNBS are obtained by using the single-mode Galerkin discretization and IHB method, which are then validated by using the numerical integration method. Moreover, Floquet theory is employed to determine the stability of obtained periodic solutions for different configurations of the nonlinear MNBS. Using the IHB method, we obtain an incremental relationship with the frequency and amplitude of time-varying axial load, which defines stability boundaries. Numerical examples show the effects of different physical and material parameters such as the nonlocal parameter, stiffness of viscoelastic medium and number of nanobeams on Floquet multipliers, instability regions and nonlinear amplitude–frequency response curves of MNBS. The presented results can be useful as a first step in the study and design of complex micro/nanoelectromechanical systems.


Introduction
Most dynamic stability studies in the literature are considering single or double-micro/nanostructure-based systems [1][2][3][4]. Investigation of the nonlinear dynamic behaviour of such systems has been an exciting perspective in the last decade due to possible applications in design procedures of different micro/nanoengineering systems [5], especially when these systems are aimed to be exploited as vibrating nanodevices such as resonators, nanosensors, or other nanoelectromechanical systems. However, it has been shown that organized nanostructure architectures made of vertically aligned forests, yarns, and sheets of carbon nanotubes (CNTs) give an exciting perspective to scale up the properties of individual CNTs and realize new functionalities [6]. The lack of reliable nonlinear dynamic models of multiple-nanostructure-based systems makes a future investigation in this field as an attractive task for researchers. A seminal idea of this work is to fill this gap in the literature and propose new models and procedures of a solution to investigate the dynamic stability of the geometrically nonlinear multiple-nanobeam system (MNBS). The presented model of MNBS consists of multiple individual simply supported nanobeams that are parallel to each other and placed into a viscoelastic medium. In general, such model can represent some nanocomposite material composed of vertically aligned CNTs array placed in some polymer. Therefore, the presented model of the nonlinear MNBS might be important to investigate the aligned arrays of CNTs [7] and CNT/polymer composites [8], especially stability of CNTs embedded in the polymer matrix. In addition, due to their remarkable physical and chemical characteristics, various applications and mathematical models of CNTs are suggested in the literature [9]. Based on experimental results, it is shown that small-scale effects play a significant role in the static and dynamic behaviour of micro/nanostructures. Since models based on classical continuum theory are scale-free, they need to be modified introducing some new assumptions in order to take into account size effects that are present at small scales. Eringen [10,11] extended the classical continuum theory by introducing integral and differential forms of constitutive equations with a single material parameter, such that it takes into account forces between atoms and internal length scale. Moreover, it has been shown that the results obtained by using the nonlocal continuum models are in good agreement with the results obtained via molecular dynamics simulations. Many authors have analysed the mechanical behaviour of different types of nanostructures within the framework of nonlocal continuum mechanics.
A pioneering study where nonlocal elasticity theory is proposed to model micro/nanostructures is the work by Pedinson et al. [12]. They analysed the static and dynamic behaviour of a simple micro/nanobeam and proposed a possible application to real MEMS and NEMS. Reddy [13][14][15] has derived new equations of motion for Euler-Bernoulli, Timoshenko, Reddy and Levinson beam theories based on the Eringen's nonlocal elasticity theory and then obtained analytical solutions for bending, vibrations, and buckling response of beams for simply supported boundary conditions. The new shear deformation beam theory proposed by Huu-Tai [16] is used to analyse the free vibration and stability of short nanobeams using Eringen's nonlocal elasticity theory. Aydogdu [17] employed the Euler-Bernoulli, Timoshenko, Levinson, Reddy, and Aydogdu beam theories to analyse the bending, buckling, and vibration of nanobeams in an analytical manner. Ansari et al. [18,19] compared results obtained from nonlocal continuum mechanics with the results obtained by molecular dynamic simulations for simple models of nanostructures and concluded that the results obtained from both theories are in good agreement. More complex nanoscale systems consist of two and more organized nanostructures such as rods, beams, and plates, which are usually coupled through some medium. The dynamic behaviour of such systems is interesting to observe and it is still not fully explored in the literature. A doublenanostructure-based system is the simplest model of a coupled multiple-nanostructure systems, which can be composed of two nanobeams, nanorods or nanoplates coupled through a medium with elastic or viscoelastic properties. Murmu and Adhikari conducted detailed vibration studies of a double-nanorod, nanobeam, and nanoplate systems [20][21][22][23], where partial differential equations of motion are obtained based on the D' Alembert's principle and nonlocal elasticity theory and solved by using analytical methods. They investigated the influence of small-scale effects and other physical parameters on natural frequencies and critical buckling loads and compared analytical results with results obtained by molecular dynamics simulations. One of the first dynamic behaviour studies of multi-layered graphene sheet system is the paper by He et al. [24] and Liew et al. [25], where the authors derived an explicit formula for natural frequencies by utilizing the classical elasticity model. Using the methodology that was proposed for chain-like mechanical systems by Rašković [26] and multiple coupled structural elements by Hedrih [27], Karličić et al. [28] presented a straightforward method to obtain analytical solutions for natural frequencies and critical buckling loads of multiple nanorods, nanobeams, and nanoplates systems based on the Eringen's nonlocal elasticity theory and trigonometric method. To this time, the mechanical behaviour of different nanostructures such as the longitudinal vibration of nanorods [29], transverse vibration of nanobeams [30][31][32], nanoplates [33] and nanoshells [34] with elastic properties are given in the literature.
Forced and parametric vibration problems of complex coupled nanostructure-based systems such as coupled nanobeams, nanorods, nanoplates have attracted much attention of the scientific community. Simsek [35] investigated the influence of a moving load on different types of single-walled carbon nanotubes. Karaoglu and Aydogdu [36] investigated the forced vibration of carbon nanotubes via nonlocal Euler-Bernoulli beam model. Ansari et al. [37,38] analysed the forced vibration of nanobeams as mechanical models of carbon nanotubes by using different numerical technics and compared the results with those obtained by the molecular dynamic simulations. In addition, Kiani [39][40][41][42] has examined the influence of different physical fields and shear effects on the free and forced vibration response of nanoscale systems composed of multiple nanobeams coupled in membrane and forest configurations. Also, Kiani [43][44][45] conducted several studies on the influence of moving loads on double-carbon nanotube system based on the nonlocal elasticity theory. Arani et al. [46] investigated the nonlinear dynamic stability of a double-graphene sheet with integrated actuators and sensors based on the Gurtin-Murdoch elasticity theory. Further, in series of papers Wang et al. [47,48] investigated the nonlinear behaviours of double-nanoplate systems for forced and parametric excitations. Recently, Pavlović et al. [49] observed the stochastic stability problem by determining the stability boundaries and regions of almost sure stability of a linear multi-nanobeam system embedded in a viscoelastic medium by using the Monte Carlo simulation method.
In general, dynamic stability analysis of nanobeamlike structures can play a significant role in design procedures of future nanodevices. For axially loaded nanobeams, where loads are time-dependent harmonic functions, a failure due to dynamic instability might occur for much smaller amplitudes of load than the failure induced by static buckling. These instability conditions usually lead to the failure of micro/nanodevices. Based on that fact, authors usually intend to analyse stability regions caused by primary parametric resonance, where the frequency of excitation is two times larger than the first natural frequencies of MNBS [50, p. 23]. Therefore, the main aim of this paper is to extend the previous studies by investigating the dynamic stability of more complex systems such as the nonlinear MNBS and to explore its primary resonance state. In addition, investigation of the influence of small-scale parameter on the stability of periodic solutions and instability regions is a very important task for the design of different types of micro/nanodevices.
It is well known that CNT/polymer composites can be observed as materials with viscoelastic properties [51][52][53]. Based on experiments and complex modulus characterization results [52], one can obtain unique properties of composite materials with aligned carbon nanotubes embedded in the polymer matrix. Knowing storage and loss modules is the first step in the identification of parameters of chosen viscoelasticity model to represent such materials. In addition, it is confirmed that in CNT/polymer composites, relaxation time and viscosity of polymers are temperature-dependent quantities [53], which can be easily obtained from complex modulus under the assumption of isothermal deformation [54]. On the other side, a model of MNBS embedded in the elastic medium is suitable to describe the van der Waals interactions between adjacent nanobeams, where values of the stiffness coefficient of the elastic medium can be determined from van der Waals forces using the methodology from the literature [39][40][41][42].
Up to this point, no investigation of the dynamic stability of the nonlinear MNBS has been carried out within the framework of nonlocal elasticity theory. In the present paper, we carry out the semi-analytical procedure based on the incremental harmonic balance (IHB) method to find periodic solutions and instability regions of the axially loaded system of m coupled nonlinear nanobeams embedded in the viscoelastic medium. We assume that nanobeams are having the same material and geometric properties as well as boundary conditions. By considering the Euler-Bernoulli beam theory, nonlocal constitutive relation and von Karman nonlinear strains, we obtain a system of m nonlinear partial differential equations of motion. Single-mode Galerkin discretization will be employed to obtain a system of m nonlinear differential equations that will be solved using the IHB method in order to obtain semi-analytical periodic solutions of the nonlinear MNBS for different configurations. Moreover, the stability of periodic solutions will be examined by introducing the Floquet theory. The parametric study will be conducted to study the influence of different parameters on the regions of instability, nonlinear amplitude-frequency response curve and Floquet multipliers. This study can be a starting point for  Let us consider the system formed by the set of straight and parallel nanobeams with geometric nonlinearity, which are embedded in a viscoelastic medium and subjected to time-dependent axial compressive loads F 1 (t) = F 2 (t) = · · · = F m (t) = F (t), Fig. 1. All nanobeams in MNBS are having the same mate-rial and geometric properties such as elastic modulus E, mass density ρ, uniform cross sections of area A and moments of inertia I . Nanobeams in MNBS are referred to as nanobeam 1, nanobeam 2 and so on until the m-th nanobeam. The transverse displacement of the i-th nanobeam is w i (x, t) , i = 1, 2, 3 . . . m. This study is limited to the case of the Euler-Bernoulli's beam model with simply supported boundary conditions (see Fig. 1). It considers only MNBS coupled in the clamped chain system, where the first and the last nanobeam in the system are coupled with a fixed base through a viscoelastic medium of stiffness k and viscosity b. Other nanobeams in MNBS are also coupled through the viscoelastic medium of stiffness per length k and viscosity b.
In order to explore the dynamic stability of nanoscale structures, we first need to consider size effects that might significantly change dynamic properties of nonlinear dynamic systems. For this purposes, modified continuum theories shown to be a reliable tool for modelling of nanoscale systems [13][14][15]. Eringen [10,11] derived a new constitutive relation in the integral form and later on in the differential form, based on the assumption that the stress at some point of an elastic body is a function of strains at all other points of that body. The aforementioned differential form of the nonlocal constitutive relation for one-dimensional case is given as where E is Young's modulus of the nanobeam, μ = (e 0 a) 2 is the nonlocal parameter (length scale) with a denoting the internal characteristic length and e 0 denoting the dimensionless material constant. Further, σ x x is the normal nonlocal stress, whereas the von Karman nonlinear strain-displacements relation is given as follows: In the literature, it is shown that the exact value of the nonlocal parameter is scattered and depends on applied boundary conditions, material and geometric properties of nanostructures. The most common method for determination of the nonlocal parameter for simple nanostructures is the fitting procedure based on molecular dynamics simulations. As mentioned previously, the nonlocal parameter μ = (e 0 a) 2 is determined by the internal characteristic length a (lattice parameter, granular size or distance between C-C bounds that is for SWCNT a =1.42 (Å), e.g. see [55]) and constant e 0 , whose value is different for each material. Eringen [10] obtained the value of nonlocal parameter by comparing the dispersion curves of plane waves from nonlocal model with those obtained by the atomic Born-Karman model of lattice dynamics. The author has found that the maximum difference of 6% is obtained for the value of constant e 0 = 0.39. Furthermore, the value of key parameter e 0 can be found by comparing the results for frequen-cies found by the nonlocal elasticity model with those from the MD simulations. In order to analyse the multiwall carbon nanotube-based system, Sudak [56] proposed that e 0 = 112.7, for the case when the, nonlocal parameter is (e 0 a) = 1.50 × 10 −8 (cm) and a = 1.42 (Å). In the paper proposed by Zhang et al. [57], the authors have estimated the value of material constant e 0 ≈ 0.82 by matching the results from nonlocal continuum model with the results from molecular dynamics simulations given in Sears and Batra [58]. Ansari et al. [19,59,60] published a series of papers where values of nonlocal parameter are calibrated based on molecular dynamic simulations. The authors have employed the NanoHive simulator to perform quasi-static molecular dynamics simulations on SWCNTs with different chirality, aspect ratios, and boundary conditions. By using the nonlinear least-square fitting procedure, fundamental frequencies are obtained for different nonlocal beam theories and matched with those calculated from MD simulations, where (e 0 a) is set as the optimization variable for each nonlocal beam model. They obtained the values of nonlocal parameter (e 0 a) from the fitting procedure for both, armchair and zigzag SWCNTs, taking into account different boundary conditions. It is shown that values of nonlocal parameter are ranging from 0.13 to 1.42 for different nanobeam theories and boundary conditions. Based on previous results, in this study we adopted that the values of nonlocal parameter goes in the range (e 0 a) = 0 − 2 (nm).
The equations of motion can be expressed in terms of displacements w i (x, t), where by taking into account nonlocal constitutive relation (1) and von Karman nonlinear strain-displacements relation (2) in the D' Alembert's principle, we obtain the following equation of motion for the i-th nanobeam [59], and clamped chain conditions as follows: where the axial load is F = F 0 + F 1 cos 2 t with amplitudes of static F 0 and dynamic part F 1 of axial load, while is the frequency of the dynamic part of axial load. From Eqs. (3) and (4), we obtain equations of motion of the clamped chain system, shown in Fig. 1, as: Mathematical expressions for initial and boundary conditions of simply supported Euler-Bernoulli nanobeams, of the same length L, are given as follows: *initial conditions

Galerkin discretization method
Now, we can reduce the system of m nonlinear partial differential equations of motion (5) to the system of m nonlinear ordinary differential equations by using the Galerkin discretization, where obtained equations represent only the time-varying part of the solution.
The first mode of approximation is assumed to be in the following form: in which r = I A is the radius of gyration of the cross section, f i (t) is the corresponding time function and φ (x) is the linear mode shape function determined from the boundary conditions. The mode shape function for corresponding boundary conditions of MNBS is given as follows: Inserting Eq. (7) into (5), multiplying the results by corresponding linear mode shape function Eq. (8), and then integrating them over the length of the nanobeam, we obtain a set of m nonlinear ordinary differential equations expressed in terms of corresponding time functions f i (t) as follows: where ω L is the natural frequency of corresponding linear system, υ is the damping ratio and γ is the nonlinear stiffness and ω 2 K , ω 2 F and are reduced constants given as follows: In the following, we introduce the IHB method to determine the periodic solutions of the presented system of m nonlinear ordinary differential equations, instability regions of the nonlinear MNBS and Floquet multipliers. It should be noted that stability of determined periodic solutions of the nonlinear MNBS is studied by using the Floquet stability theory.

Instability regions and periodic solution: IHB method
Lau et al. [61] developed the methodology for the IHB method which is later on applied to different problems [62][63][64][65][66]. In general, IHB is the semi-analytical method that can be used for finding the solutions of differential equations appearing in mathematical models of various problems in mechanics. In order to obtain instability regions and periodic solutions of the system of nonlinear ordinary differential equations (9), we use the IHB method to obtain iterative relationships with frequency, response amplitudes and excitation amplitude of time-varying axial load. First, one should introduce a new time scale τ = t into Eq. (9) to obtain the system of nonlinear ordinary differential equations in the following form: Dividing Eq. (15) with ω 2 L , one can obtain a new form of Eq. (11) as follows: where the values of the frequency depends on the configuration of the nonlinear MNBS and˜ = and ˜ are corresponding increments of amplitude, frequency and amplitude of excitation load, respectively, for neighbouring vibration state, which is given as follows: Substituting relations (13) into Eq. (12) and neglecting higher-order terms, a linearized incremental relation can be expressed as Assuming the periodic solutions of Eq. (14), the functions f i0 (τ ) and f i (τ ) can be taken as Fourier series in the following forms: where Inserting Eq. (15) into the system of Eq. (14) and applying the Galerkin procedure presented in [61][62][63][64][65][66][67][68], we obtain the system of linear algebraic equations in terms of unknown increments of amplitudes A i in the following form: . . . where In Eq. (17), corrective vector term [R] {A} tends to zero when the values of f i0 (τ ) , 0 and˜ 0 tends to exact solutions. Set of linearized algebraic equations presented in Eq. (17) can be solved incrementally by starting from any known linear solution using -incrementation procedure.
Since the system of Eq. (17) has two more unknowns than the number of equations in the system, we need to set one value of Fourier coefficients in Eq. (15) to be constant and corresponding increment equal to zero. In the˜ -incrementation process, the ˜ is an active increment whereby giving the initial values of A i , 0 and˜ 0 one can obtain principal instability regions that are bounded by two boundaries. When using the Newton-Raphson method for solving the system (17), we need to do the following two steps. First, we determine the initial values of A i , 0 and 0 , set a i1 = Const., b i1 = Const., a i1 = 0, b i1 = 0 and consider ˜ as an active increment during simulation to obtain the same number of equations and unknowns in Eq. (17). Second, by using the iterative Newton-Raphson method, we solve the incremental relation (17), where the results for A i and are obtained iteratively until the residue Euclidian norm | [R] {A} | is smaller than a preset tolerance for In the next step, the initial value of˜incrementation procedure is˜ 0 + ˜ , and we repeat an iterative procedure until the residue Euclidian norm | [R] {A} | is smaller then a preset tolerance. The process is repeated until˜ 0 reaches the set value.
In order to obtain periodic solutions of the system of nonlinear ordinary differential equations (12), we reduce the incremental relations to and for a periodic solution we obtain where Then, the incremental relation (17) becomes

Stability of the periodic solution: Floquet theory
In order to analyse the stability of periodic solutions of the system of nonlinear ordinary differential equations (9), we consider the Floquet theory in R N . Now, we write the system of nonlinear ordinary differential equations for nonlinear MNBS in general form as follows: where Introducing small perturbations f (τ ) in a neighbourhood of the periodic solution f 0 (τ ), i.e. by letting we can analyse the stability of periodic solutions by analysing the system of linear equations with variable coefficients in terms of small perturbations f (τ ). Introducing Eq. (22) into the system of Eq. (21) and after linearization, we obtain the system of linear differential equations as follows: (23) where f 0 (τ ) is the periodic solution determined by considering the IHB method. It should be noted that Eq. (23) represents the system of perturbed equations around known solutions f 0 (τ ). The stability characteristics of known periodic solutions are found through the multi-variable Floquet theory.
We introduce state-space variables where for Y (τ ) = f , f T and G (τ ) denoting the periodic matrix. Based on the Floquet theory [68][69][70], the stability criteria for determined periodic solutions of MNBS is related to eigenvalues of the matrix M (Floquet multipliers), which is the transition matrix and the real parts of characteristic exponents (Floquet exponents). In the case when all values of Floquet multipliers are located inside the unit circle centered at the origin of the complex plane, the periodic solutions are stable or asymptotically stable. However, when the values of Floquet multipliers are out of the unit circle, then the periodic solutions are unstable [68][69][70].
In order to numerically determine the Floquet multipliers, we assume that period T = 2π of f 0 (τ ) is divided into N k subintervals, in which the k-th interval is k = τ k − τ k−1 for τ k = kT /N k . Moreover, G (τ ) is the continuous matrix with respect to τ , so in the k-th interval it can be replaced by the constant matrix provided in the case when N k is chosen to be sufficiently large Now, we can write the transition matrix in the following form: where N j denotes the number of terms in the approximation of the constant matrix G k . From the transition matrix M, one can obtain Floquet multipliers as eigenvalues of Eq. (26) in the following form: For the nonlinear MNBS consisting of m nanobeams, the periodic matrix G (τ ) has the following form: where P i , (i = 1, 2, . . . m), R and T are defined as {a ik cos (kτ ) From the relation (28), it is easy to obtain a system matrix for a specific configuration of the nonlinear MNBS.

Numerical results
This section is divided into four parts. First, we show a comparative study where results for periodic solutions of the system of nonlinear differential equations obtained by employing the IHB method are validated The system of m nonlinear differential equations (9) represents a mathematical model of the nonlinear MNBS that is vibrating under the influence of timevarying axial loads with considered initial conditions. In the case of larger systems, such as the nonlinear MNBS, it is difficult to find system responses using only analytical methods. In addition, in a large number of nanoengineering applications, the amplitude of vibration is not small and for such systems, we can only use numerical or approximate technics. It should be noted that IHB method has several advantages in comparison with classical approximated methods. First, there are no limitations such as small exciting parameters or weak nonlinearities and second, it is easy for computer implementation. Because of these facts, we determine periodic solutions of the nonlinear MNBS without introducing any limitations such as the bookkeeping parameter.

Periodic solutions and stability
In this study, we consider several different cases of MNBS with m = 1, 3, 5 nanobeams and for different initial conditions. We adopted the material and geometric properties for nanobeams in MNBS corresponding to the properties of single-walled carbon nanotubes as: Young's modulus E = 1.1 TPa, density ρ = 1300 kg m 3 , length L = 45 nm, inner d 1 = 3 nm and outer d 0 = 2.32 nm diameters of the nanotube. For the numerical purposes, we use the value of elastic coefficient of the medium as k = 10 8 (N/m 2 ), while influence of damping coefficient is almost neglected and small value of b = 10 −6 (Ns/m 2 ) is adopted.
In order to illustrate the accuracy of the proposed IHB methodology, we compare the obtained results with the results from the Runge-Kutta method in Figs. 2, 3 and 4, where a fine agreement of the results is achieved. These figures show time response functions of the first nanobeam and periodic orbits of each nanobeam for different configurations of the nonlinear MNBS. Figure 2 shows periodic orbits of the special case of the nonlinear MNBS consisting of a single nanobeam, which is the case very often analysed in the literature. However, in Figs. 3 and 4, one can observe periodic orbits of the nonlinear MNBS consisting of three-and five-coupled nanobeams. We can notice that Now, we can investigate the stability of determined periodic solutions of the nonlinear MNBS subjected to axial time-dependent loads (Figs. 5, 6, 7). By considering the Floquet theory, we can obtain Floquet multipliers for different configurations of the nonlinear MNBS consisting of one, three and five nanobeams in the system. By calculating the Floquet multipliers and tracing their evaluation, we predict the stability or identify bifurcations of the presented system. For stability of the periodic solution, Floquet multipliers must be centred within the unit circle with the origin in the centre of the complex plane. However, if the values of Floquet multipliers leave the unit circle, we can predict bifurcation according to [71]. For particular values of material parameters and initial conditions, we consider the iterative equation (20) and obtain the following stable periodic solutions wherein the stability of the trajectory is confirmed by using the Floquet theory and above-described iterative procedure [69,70]. We set the values of N k = 10 4 and N j = 5 in the relation (26) showing Floquet multipliers as functions of amplitudes of axial loads and nonlocal parameter, Figs. 5, 6 and 7, for different configurations of the nonlinear MNBS.
The influence of axial load amplitude on the Floquet multipliers for the special case of MNBS consisting of a single nanobeam is shown in Fig. 5. Moreover, we get Floquet multipliers as complex numbers, where real parts are shown in Fig. 5a and imaginary parts in Fig. 5b.   By increasing the value of axial load˜ , both parts of the Floquet multipliers remain within the unit circle of the complex plane, as shown in Fig. 5c. Based on that fact, we can conclude that obtained periodic solutions for the nonlinear system are stable for both values of the nonlocal parameter. It should be noted that the real parts of Floquet multipliers decrease for an increase in the nonlocal parameter, while the imaginary parts increases. Figure 6 shows three values of Floquet multipliers as functions of the nonlocal parameter and amplitude of axial load for the nonlinear MNBS consisting of three nanobeams. From obtained curves, it can be observed that the influence of the amplitude of axial load on the Floquet multipliers is not a smooth function coming out of the unit circle in the complex plane, as shown in Fig. 6g, i. Based on that fact, we can conclude that the periodic solution of the third nanobeam in MNBS is unstable. For the Floquet multipliers regarding the first and the second nanobeam, we can observe that obtained periodic solutions are stable, as shown in Fig. 6c, f. It should be noted that the nonlocal parameter has a significant influence on the first two Floquet multipliers while for the third one it almost vanishes, as shown in Fig. 6i.
The Floquet multipliers determining the stability of the nonlinear MNBS composed of five nanobeams are shown in Fig. 7. We can observe a very interesting behaviour of Floquet multipliers with a significant dependence on the nonlocal parameter. From obtained results, we can notice that only the fourth and the fifth nanobeam in MNBS are having stable periodic solutions with the values of Floquet multipliers remaining within the unit circle of the complex plane, as shown in Fig. 7l, p. However, other periodic solutions of MNBS are unstable, where the values of the Floquet multipliers are out of the unit circle, as shown in Fig. 7c, f, i. From the physical point of view, we can conclude that the number of nanobeams decreases the overall stiffness of the MNBS and this fact lead to the reduction of the system stability. where the value of the increment is λ = 0.001, with five cosine and five sine harmonic terms. The boundary lines leaning to the left are determined by using the sine terms, while boundary lines leaning to the right are defined by the cosine terms only. When considering incremental relations to determine instability regions, we can examine cases with small and large initial amplitudes as well. Small initial amplitudes lead to the linear case while large values of initial amplitudes lead to the nonlinear case, for the same values of the excitation force amplitude. Both cases are explored, where amplitudes are shown as functions of the amplitude of axial load, a number of nanobeams and different physical parameters such as stiffness coefficients (Fig. 8), the nonlocal parameter (Fig. 9) and the amplitude of static axial load (Fig. 10). The values of initial amplitudes for the incremental equation (17) are given as a i1 = 0.1 and b i1 = 0.1 for the linear case and a i1 = 3 and b i1 = 3 for the nonlinear case. Applying the Newton-Raphson procedure, where incrementation process starts from λ = 0 to λ = 1, we can investigate the instability region of the nonlinear MNBS. For the present analysis, we adopted the same material and geometrical properties of nanobeams as in the previous case. The instability regions are represented by the areas that are surrounded by two stability boundaries. Figure 8 shows the influence of stiffness coefficient k and the number of nanobeams on the instability regions of the linear and nonlinear cases. It can be concluded that an increase in stiffness coefficient k leads to narrowing of instability regions while an increase in the number of nanobeams in MNBS leads to an increase in instability regions. From the physical point of view, an increase in the number of nanobeams decreases the overall stiffness of the system and dynamic system is becoming more unstable. From the nonlinear case (Fig. 8b), we can notice that different configurations of MNBS shifts the instability regions, which is the consequence of the system's nonlinearity.
The influence of the nonlocal parameter on instability regions of MNBS for different configurations is shown in Fig. 9. One can observe that the influence of the nonlocal parameter on instability regions is small for small values of the amplitude of axial load in the linear case. However, in the nonlinear case, the nonlocal parameter has a significant influence on instability regions. As expected, we have shifts of instability regions of MNBS for a different number of nanobeams in the system. The amplitudes of nonlinearity and over-  The last Fig. 10 illustrates the effects of the amplitude of static axial load on instability regions of MNBS. An increase in the amplitude of static axial load increases the instability regions. This increase is smaller for the nonlinear case than for the linear one. In addition, this effect occurs as a result of an increase in the overall stiffness of the system. For all the presented cases, an increase in the amplitude of axial load˜ leads to an increase in instability regions and the system becomes more unstable.
In a general case, we can conclude that change of a number of nanobeams in MNBS affects its instability regions in such a manner that these regions are expanding for an increase in the parameter m and small initial amplitudes. On the other side, for larger initial amplitudes, we can notice that the effect of a number of nanobeams goes in two directions. First, the instability regions are shifted relative to the case of small initial amplitudes and second, instability regions are increasing for a change of a number of nanobeams in the system. From the mechanical point of view, we can conclude that the overall stiffness of MNBS is reduced when increasing the number of nanobeams, which consequently results in reduced stability of the system.

The free nonlinear vibration of MNBS
In order to analyse the free nonlinear vibration of MNBS, amplitude-frequency responses for different configurations of MNBS obtained by using the IHB method are presented in Figs. 11 and 12. It should be noted that nonlinear amplitude-frequency curves are determined by using the b 1 -incrementation, where the value of the increment is b 1 = 0.001, using the five odd cosine terms [67]. In this subsection, we will neglect the influence of the axial load and damping coefficient in Eq. (5) to obtain the governing equation for the free vibration of the nonlinear MNBS embedded in the elastic medium. Figure 11 shows the influence of the stiffness coefficient of the elastic medium on the nonlinear amplitude-frequency response curve of MNBS consisting of one-, three-and five-coupled nanobeams. It can be observed that an increase in the stiffness coefficient tends to decrease nonlinear effects visible on amplitude-frequency response curves, while these effects almost vanish for the highest values of stiffness coefficient. Moreover, in the case of the larger stiffness of the elastic medium, nonlinear frequencies are approaching the linear ones. On the other hand, an increase in the number of nanobeams in MNBS leads to a decrease in overall stiffness and nonlinearity becomes more prominent.
The effects of the nonlocal parameter and a number of nanobeams on the nonlinear amplitude-frequency response curves for each case of nanobeams in the system are shown in Fig. 12. It can be noticed that for given amplitudes and for an increase in the nonlocal parameter, the frequency ratio increases but the frequency decreases. This fact leads to the conclusion that the nonlocal parameter softens the system, i.e. reduces the overall stiffness. Moreover, one can observe that the influence of nonlocal effects is larger for a higher number of nanobeams and larger vibration amplitudes.

Validation study
In this work, results are obtained for MNBS embedded in the viscoelastic medium and the authors did not find a similar mechanical model in the literature that can be used for the validation study. However, for the special case of MNBS, i.e. the case with a single nanobeam on elastic foundation, we use the paper by Fu et al. [67]. Figure 13 shows the nonlinear amplitude-frequency response of MNBS consisting of one, three and five nanobeams and the amplitude-frequency response of  [67] for (e 0 a) = 1 nm the system with a single nanobeam resting on the elastic foundation presented in Fu et al. [67]. It should be noted that our model of MNBS is coupled in clamped chain system, i.e. for the case of a single nanobeam it is coupled with the fixed base through the elastic medium on both sides. This leads to higher values of the overall stiffness and such system is more constrained than the model proposed by Fu et al. [67].

Scope and limitations of the presented model
The given nonlinear MNBS is based on the nonlocal Euler-Bernoulli beam model and von Karman nonlinear deformations and does not take into account shear deformations. However, in order to consider this effect, one should consider Timoshenko or higher-order shear deformation beam theories. From the literature [18], we can see that the effect of transverse shear deformation leads to lower vibration frequencies comparing to the case with Euler-Bernoulli beam. This effect is more prominent in higher vibration modes. Therefore, the consideration of transverse shear deformation leads to a decrease in stability of the presented system, i.e. it leads to widening of instability regions [72]. It should be noted that applied axial loads are smaller than the critical buckling load, and MNBS is not buckled. However, the nonlinear vibration behaviour of buckled MNBS is also a very important issue and should be studied in some future works.
The second limitation refers to the fact that presented MNBS model consider amplitudes of nanobeams that are bounded by the surrounding viscoelastic medium. Based on that, we analysed only small nonlinear deformations by using the von Karman theory. It should be emphasized that we employed the phenomenological viscoelastic model based on simple mechanical analogy to describe the rheological behaviour of the viscoelastic medium in MNBS. Since the applied forcedisplacement relation is a phenomenological model of Kelvin-Voight's viscoelastic stress-strain constitutive relation and reaction force of viscoelastic medium is distributed over nanobeam's length, it can be represented by parallelly connected springs and dashpots with corresponding elasticity and viscosity coefficients. An interested reader can find more on experimental works and dynamic behaviour of CNT/polymer nanocomposites in [51][52][53]75].
In the literature, the authors have found a series of papers [61,67,73,74] where linear mode shape functions are used for a nonlinear vibration analysis of classical and scale-dependent beam models. However, we should mention two interesting papers by El-Borgi et al. [73] and Nayfeh and Lacarbonara [74], where the authors concluded that the Galerkin single-mode method is suitable for the discretization of governing equations of nonlinear vibration of simply supported beams, especially for the primary resonance state and the first vibration mode. According to this, the authors of the present paper employed only a single-mode Galerkin method for discretization of the governing equations of given MNBS. One should know that each nanobeam in MNBS has the same simply supported boundary conditions and therefore, same linear mode shape functions are considered in the assumed solutions of the system of governing equations. Moreover, in the absence of internal resonances, the steady-state response contains only the directly excited modes. The mathematical model of MNBS with considered interactions between modes, internal and combined resonances, is interesting to be analysed in some future study.

Conclusions
The present analysis describes the dynamic stability problem of the nonlinear MNBS embedded in a viscoelastic medium, where each nanobeam in the system is simply supported and subjected to axial time-dependent loads. A mathematical model of the nonlinear MNBS is based on the system of nonlinear partial differential equations of motion derived by using the Eringen's nonlocal elasticity theory, Euler-Bernoulli's beam theory and nonlinear von Karman strain-displacement relation. By employing the Galerkin discretization technique and IHB method, we determined the stability boundaries and periodic solutions of the system of nonlinear ordinary differential equations. For such determined periodic solutions, we introduced a Floquet theory to analyse their stability. In the numerical study, we have shown the influence of different material and geometric parameters on the stability behaviour of the presented nonlinear MNBS. A few novelties are proposed in this paper: • Semi-analytical periodic solutions are obtained for nonlinear MNBS by using the IHB method, without introducing small parameter.
It can be concluded that the analysis of periodic solutions of the nonlinear MNBS without taking into account nonlocal effects might lead to large errors when investigating the dynamic behaviour of such complex systems. It is well known that a very small variation in amplitudes of axial loads can cause structural collapse and the values of amplitudes should be chosen carefully. Information based on reliable models are important for the design of future nanostructurebased devices; therefore, this study besides theoretical could also have practical applications in nanoengineering fields. In addition, this study might be important for future theoretical analyses of the nonlinear dynamic behaviour of MNBS such as internal resonance, combined resonance conditions, and chaos.