Introducing dimensionless variables in Eq. (19) and Eq. (33) is convenient. Lets \(\bar {x}=x/l\), \(\bar {w}=w/l\), \(\bar {x}=x/l\), \(\bar {t}=t/{\omega _0}\), and \(\omega _{0}^{2}={\pi ^4}EI{\left( {{l^4}m} \right)^{ - 1}}\). Assumes the loads are uniform, namely \(\bar {F}\left( x \right)=\text{c}\text{o}\text{n}\text{s}\text{t}\), so \({\partial ^2}\bar {F}/\partial {x^2}=0\). Eq. (19) for the extensibility effect can rewrite as

$$\begin{gathered} \frac{{{\partial ^2}}}{{\partial {{\bar {t}}^2}}}\left( {\bar {w} - \frac{{{\mu ^2}}}{{{l^2}}}\frac{{{\partial ^2}\bar {w}}}{{\partial {{\bar {x}}^2}}}} \right)+\frac{C}{{m{\omega _0}}}\frac{{\partial \bar {w}}}{{\partial t}}+\frac{1}{{{\pi ^4}}}\frac{{{\partial ^4}\bar {w}}}{{\partial {{\bar {x}}^4}}} \hfill \\ +\frac{{{N_0}}}{{m\omega _{0}^{2}{l^2}}}\left( {\frac{{{\partial ^2}\bar {w}}}{{\partial {{\bar {x}}^2}}}\left. { - \frac{{{\mu ^2}}}{{{l^2}}}\frac{{{\partial ^4}\bar {w}}}{{\partial {{\bar {x}}^4}}}} \right) - \frac{{2\lambda ID}}{{m\omega _{0}^{2}{l^6}}}\frac{{{\partial ^2}}}{{\partial {{\bar {x}}^2}}}{{\left( {\frac{{{\partial ^2}\bar {w}}}{{\partial {{\bar {x}}^2}}}} \right)}^3}} \right. \hfill \\ - \frac{{EA}}{{2m\omega _{0}^{2}{l^2}}}\left( {\frac{{{\partial ^2}\bar {w}}}{{\partial {{\bar {x}}^2}}}} \right.\left. { - \frac{{2\lambda +{\mu ^2}}}{{{l^2}}}\frac{{{\partial ^4}\bar {w}}}{{\partial {{\bar {x}}^4}}}} \right)\int\limits_{0}^{1} {\left[ {{{\left( {\frac{{\partial \bar {w}}}{{\partial \bar {x}}}} \right)}^2}+\frac{{2\lambda }}{{{l^2}}}{{\left( {\frac{{{\partial ^2}\bar {w}}}{{\partial {{\bar {x}}^2}}}} \right)}^2}} \right]d{\kern 1pt} \bar {s}} =\frac{{\bar {F}}}{{m\omega _{0}^{2}l}}. \hfill \\ \end{gathered}$$

35

Eq. (33) for inextensible beams can rewrite as

$$\begin{gathered} \frac{{{\partial ^2}\bar {w}}}{{\partial {{\bar {t}}^2}}} - \frac{{{\mu ^2}}}{{{l^2}}}\frac{{{\partial ^4}\bar {w}}}{{\partial {{\bar {t}}^2}\partial {{\bar {x}}^2}}}+\frac{C}{{m{\omega _0}}}\frac{{\partial \bar {w}}}{{\partial t}}+\frac{1}{{{\pi ^4}}}\frac{{{\partial ^4}\bar {w}}}{{\partial {{\bar {x}}^4}}}+\frac{{{N_0}}}{{m\omega _{0}^{2}{l^2}}}\left( {\frac{{{\partial ^2}\bar {w}}}{{\partial {{\bar {x}}^2}}} - \frac{{{\mu ^2}}}{{{l^2}}}\frac{{{\partial ^4}\bar {w}}}{{\partial {{\bar {x}}^4}}}} \right) \hfill \\ - \frac{{{\mu ^2}{N_0}}}{{2m\omega _{0}^{2}{l^4}}}\frac{\partial }{{\partial \bar {x}}}\left[ {{{\left( {\frac{{\partial \bar {w}}}{{\partial \bar {x}}}} \right)}^2}\frac{{{\partial ^3}\bar {w}}}{{\partial {{\bar {x}}^3}}}} \right]+\frac{{EI}}{{2m\omega _{0}^{2}{l^4}}}\frac{\partial }{{\partial \bar {x}}}\left[ {\frac{{\partial \bar {w}}}{{\partial \bar {x}}}{{\left( {\frac{{{\partial ^2}\bar {w}}}{{\partial {{\bar {x}}^2}}}} \right)}^2}} \right. \hfill \\ \left. {+{{\left( {\frac{{\partial \bar {w}}}{{\partial \bar {x}}}} \right)}^2}\frac{{{\partial ^3}\bar {w}}}{{\partial {{\bar {x}}^3}}}} \right] - \frac{{DI}}{{m\omega _{0}^{2}{l^4}}}\left[ {{{\left( {\frac{{{\partial ^2}\bar {w}}}{{\partial {{\bar {x}}^2}}}} \right)}^3} - \frac{{{\mu ^2}}}{{{l^2}}}\frac{{{\partial ^2}}}{{\partial {{\bar {x}}^2}}}{{\left( {\frac{{{\partial ^2}\bar {w}}}{{\partial {{\bar {x}}^2}}}} \right)}^3}} \right]=\frac{{\bar {F}}}{{m\omega _{0}^{2}l}}. \hfill \\ \end{gathered}$$

36

For simplicity, we use hinged-hinged and simply supported beams as examples to discuss the extensible and inextensible beams correspondingly. Hence their normalized boundary conditions are identical:

$$\bar {w}\left( {0,\bar {t}} \right)=\bar {w}\left( {1,\bar {t}} \right)=\frac{{{\partial ^2}\bar {w}}}{{\partial {{\bar {x}}^2}}}\left( {0,\bar {t}} \right)=\frac{{{\partial ^2}\bar {w}}}{{\partial {{\bar {x}}^2}}}\left( {1,\bar {t}} \right)=0$$

37

.

It is challenging to solve nonlinear Eq. (35) or Eq. (36) accurately. There are two common approaches to solving a nonlinear partial differential equation approximately. One of them is that the partial differential equation is reduced to nonlinear ordinary differential equations through Galerkin method, and then the equations are solved using perturbation methods [34–37]. The second is to solve the nonlinear partial differential equation using perturbation methods directly. The direct multiscale method is an essential improvement of the classical perturbation methods, and it has been widely used to directly solve nonlinear partial differential equations more accurately than the first methods [37]. Nayfeh [38–40], Luongo [41, 42], and Lacarbonara [43–46] and their collaborators’ works may be substantial for the direct multiscale method. We apply the Galerkin and the multi-scale methods to solve the equations for simplicity. In fact, the two methods have widely applied to nonlinear equations of structures [18, 20–22, 47–50]. Under boundary conditions Eqs. (37), the approximate solutions of Eq. (35) and Eq. (36) can be written as

$$\bar {w}=\sum\limits_{{n=1}}^{\infty } {{{\bar {\eta }}_n}\left( {\bar {t}} \right)\sin n\pi \bar {x}}$$

38

.

A set of ordinary differential equations can obtain through the Galerkin truncation [34, 36]: substitutes Eq. (38) into Eq. (35) or Eq. (36) respectively, then \(\sin \left( {n\pi \bar {x}} \right)\) is multiplied by both sides of the equations, and integrates in the interval \(\left[ {0,1} \right]\). For simplicity, we only takes the first term of Eq. (38) and let \({\bar {\eta }_1}=\eta\), get

$${m_j}\ddot {\eta }+{c_j}\dot {\eta }+{k_j}\eta +{d_j}{\eta ^3}=f,\quad j=1,2.$$

39

Here \(j=1\) indicates the hinged-hinged beams, and \(j=2\) indicates the simply supported beams. The parameters in Eqs. (39) are

$$\begin{gathered} {m_1}={m_2}=1+\frac{{{\pi ^2}{\mu ^2}}}{{{l^2}}},\quad {c_1}={c_2}=\frac{C}{{m{\omega _0}}},\quad {k_1}={k_2}=1 - \frac{{{N_0}{\pi ^2}}}{{m\omega _{0}^{2}{l^2}}}\left( {1+\frac{{{\pi ^2}{\mu ^2}}}{{{l^2}}}} \right), \hfill \\ {d_1}=\frac{{{\pi ^4}EA}}{{4m\omega _{0}^{2}{l^2}}}{\left( {1+\frac{{2{\pi ^2}\lambda }}{{{l^2}}}} \right)^2} - \frac{{3{\pi ^8}\lambda ID}}{{2m\omega _{0}^{2}{l^6}}}+\frac{{{\pi ^4}EA{\mu ^2}}}{{4m\omega _{0}^{2}{l^4}}}\left( {1+\frac{{2{\pi ^2}\lambda }}{{{l^2}}}} \right), \hfill \\ {d_2}= - \frac{{3{\pi ^6}{\mu ^2}{N_0}}}{{8m\omega _{0}^{2}{l^4}}}+\frac{{{\pi ^6}EI}}{{6m\omega _{0}^{2}{l^4}}}+\frac{{3{\pi ^6}DI}}{{4m\omega _{0}^{2}{l^4}}}\left( {1+\frac{{{\pi ^2}{\mu ^2}}}{{{l^2}}}} \right),\;F=\frac{{4\bar {F}}}{{\pi m\omega _{0}^{2}l}}. \hfill \\ \end{gathered}$$

40

Eqs. (40) indicate that the linear coefficients of the two models are the same, but the nonlinear coefficients are different. Moreover, the non-local effect and the material nonlinearity are coupled in nonlinear terms. \({d_1}\) demonstrates that the influence of the material nonlinearity will decrease accompanying the increase of beam’s length, as hinged-hinged SWCNTs shown in Fig. 2. If we neglect the nonlocal effect, the material nonlinearity’s influence on \({d_2}\) is invariant accompanying the length’s change due to \(\omega _{0}^{2}{l^4}={\pi ^4}EI/m\), as simply supported SWCNTs shown in Fig. 9.

The nonlinear terms, which are induced by the finite deformations and the non-local effect, produce a hardening effect in two models due to \({\mu ^2}>0\). On the other hand, the nonlinear terms induced by the material nonlinearity have a softening effect in two models due to \(D<0\) for SWCNTs. It is interesting for hinged-hinged SWCNTs that the hardening effect will be smaller than the softening effect for the small-length tubes, as shown in Fig. 2. The following discussions will show that this competitive relationship between the softening and the hardening impacts the mechanical behavior of SWCNTs significantly.

Neglecting the inertia and damping terms in Eqs. (39), static bending deformations of the beams’ middle points are obtained as

$${k_j}\eta +{d_j}{\eta ^3}=F,\quad j=1,2.$$

41

In the following, we will research the mechanical behaviors of hinged-hinged beams and simply supported beams according to Eq. (39) and Eq. (41), respectively.

We rewrite Eqs. (39) as

$$\ddot {\eta }+{\bar {c}_j}\dot {\eta }+\bar {\omega }_{j}^{2}\eta +{\bar {d}_j}{\eta ^3}=\bar {f},\quad j=1,2$$

42

,

here \({\bar {c}_j}={c_j}/{m_j}\), \({\omega ^2}={k_1}/{m_j}\), \({\bar {d}_j}={d_j}/{m_j}\), \(\bar {f}={\bar {f}_j}=F/{m_j}\). The multiple scale method [36], which is widely used to solve weak nonlinear differential equations for macro-structures, will be applied to solve Eq. (42). Lets \({\bar {c}_j}={\varepsilon ^2}{c_j}\), \(\bar {f}={\varepsilon ^3}f\cos \left( {\Omega \bar {t}} \right)\), gives

$$\ddot {\eta }+\omega _{j}^{2}\eta +{\varepsilon ^2}{c_j}\dot {\eta }+{\bar {d}_j}{\eta ^3}={\varepsilon ^3}f\cos \left( {\Omega \bar {t}} \right),\quad j=1,2.$$

43

.

Supposes

$$\eta \left( {t;\varepsilon } \right)=\varepsilon {\eta _1}\left( {{T_0},{T_2}} \right)+{\varepsilon ^3}{\eta _3}\left( {{T_0},{T_2}} \right)$$

44

,

and substitutes it into Eq. (43), then equates the coefficients of \(\varepsilon\) and \({\varepsilon ^3}\) on both sides, has

$$\varepsilon :\quad D_{0}^{2}{\eta _1}+\omega _{i}^{2}{\eta _1}=0,$$

45

$$\begin{gathered} {\varepsilon ^3}:\quad D_{0}^{2}{\eta _3}+\omega _{j}^{2}{\eta _3}= - 2{D_0}{D_2}{\eta _1} \hfill \\ \quad - 2{c_j}{D_0}{\eta _1} - {{\bar {d}}_j}\eta _{1}^{3}+\frac{1}{2}f\exp \left( {i\Omega {T_0}} \right), \hfill \\ \end{gathered}$$

46

where \({D_0}=d/d{T_0}\), \({D_1}=d/d{T_1}\) and \(D_{0}^{2}={d^2}/dT_{0}^{2}\). The solution of Eq. (45) is

$${\eta _1}=A\left( {{T_2}} \right)\exp \left( {i{\omega _j}{T_0}} \right)+CC$$

47

,

here, *CC* means the complex conjugate. Substituting Eq. (47) into Eq. (46), gets

$$\begin{gathered} D_{0}^{2}{\eta _3}+\omega _{j}^{2}{\eta _3}=\frac{1}{2}f\exp \left( {i{\kern 1pt} \Omega {T_0}} \right) \hfill \\ - \left[ {i2\omega \left( {A^{\prime}+{c_j}A} \right)} \right.\left. {+3{{\bar {d}}_j}{A^2}\bar {A}} \right]\exp \left( {i{\omega _j}{T_0}} \right)+NST. \hfill \\ \end{gathered}$$

48

Here the prime denotes the derivative with respect to \({T_2}\), and NST denotes non-secular terms [36]. When the load’s frequency \(\Omega\) approaches the nanobeam’s modal frequency \({\omega _j}\) (primary resonance), the beam will appear a relatively large amplitude response. Under this condition, lets \(\Omega ={\omega _j}+{\varepsilon ^2}\sigma\), so the solvable condition of Eq. (48) are

$$- i2{\omega _j}\left( {A^{\prime}+cA} \right) - 3{d_j}{A^2}\bar {A}+\frac{1}{2}f\exp \left( {i\sigma {T_2}} \right)=0,\quad j=1,2.$$

49

.

Letting \(A=\left( {\alpha /2} \right)\exp \left( {i\beta } \right)\), and substituting it into Eq. (49), then separating the real part and the imaginary part, gets

$$\begin{gathered} \alpha ^{\prime}= - {c_j}\alpha +\frac{f}{{2\omega }}\sin \gamma , \hfill \\ \alpha \gamma ^{\prime}=\sigma \alpha - \frac{3}{{8{\omega _i}}}{{\bar {d}}_j}{\alpha ^3}+\frac{f}{{2{\omega _j}}}\cos \gamma , \hfill \\ \end{gathered}$$

50

here \(\gamma =\sigma {T_2} - \beta\). Steady-state motions occur when \(\alpha ^{\prime}=\gamma ^{\prime}=0\), which corresponds to the singular points of Eqs. (50). The steady-state solutions can obtain from the following algebraic equations [36]

$$\left[ {{c^2}+{{\left( {\sigma - \frac{{3{{\bar {d}}_j}}}{{8\omega _{j}^{2}}}{\alpha ^2}} \right)}^2}} \right]{\alpha ^2}=\frac{{{f^{{\kern 1pt} 2}}}}{{4\omega _{j}^{2}}},\quad j=1,2$$

51

.

The stability of the steady-state solutions is judged by investigating the nature of the singular points of Eqs. (50). Lets \(a={a_0}+{a_1}\) and \(\gamma ={\gamma _0}+{\gamma _1}\), and substitutes them into Eqs. (50), expands for small \({a_1}\) and \({\gamma _1}\), and keeps linear terms in \({a_1}\) and \({\gamma _1}\), then gets

$$\begin{gathered} {{\alpha ^{\prime}}_1}= - {c_j}{\alpha _1}+\frac{{{\gamma _1}f\cos {\gamma _0}}}{{2{\omega _j}}}, \hfill \\ {{\gamma ^{\prime}}_1}= - \left( {\frac{{f\cos {\gamma _0}}}{{2{\omega _j}\alpha _{0}^{2}}}+\frac{{3{{\bar {d}}_j}{\alpha _0}}}{4}} \right){\alpha _1} - \frac{{{\gamma _1}f\sin {\gamma _0}}}{{2{\omega _j}{\alpha _0}}}. \hfill \\ \end{gathered}$$

52

Here, it is used in Eqs. (52) that \({a_0}\) and \({\gamma _0}\) are the singular points of Eqs. (50). The stability of steady-state motions depends on the coefficient matrix’s eigenvalues of Eqs. (52). If the real parts of eigenvalues are greater than zero, the solutions are unstable [36]. Hence the steady-state motions are unstable if

$$c_{j}^{2}+\left( {\frac{{3{d_j}\alpha _{0}^{2}}}{{8{\omega _j}}} - \sigma } \right)\left( {\frac{{9{{\bar {d}}_j}\alpha _{0}^{2}}}{{8{\omega _j}}} - \sigma } \right)<0$$

53

.

Eq. (53) indicates that the nonlinear term affects the stability of the steady-state solutions. So both the nonlocal effect and the material nonlinearity affect the stability of the solutions. In the response curves in the next section (Fig. 5–7 and Fig. 11–12), the solid line represents the stable solutions and the dashed line represents the unstable solutions.