Nonlinear dynamics of a planar beam–spring system: analytical and numerical approaches

Abstract

The multiple timescales method is applied to the exact partial differential equations of the planar motion of a hinged–simply supported beam with a linear axial spring of arbitrary stiffness. The forced-damped and free oscillations of the system around frequencies corresponding to nth natural bending mode are examined thoroughly and compared with numerical simulations as well as with already published results obtained by Lindstedt–Poincaré method. A special numerical technique using explicit finite element method to draw the frequency–response curves is appositely developed. The well-known jump phenomena between resonant and non-resonant branches, as well as superharmonic resonances, have been detected numerically.

1 Introduction

Modern technologies require more performance mechanical systems, lighter and more slender than the previous one, that are subjected to large displacements and thus behave in the nonlinear regime. Thus, recent studies pay more and more attention to nonlinear dynamics of real systems. This interest is caused mainly by the need of detecting multiple solutions, for their avoidance when they represent dangerous conditions for structures, for example in bridges, towers and aeroplanes. On the other hand, nonlinear phenomena can be exploited in order to harvest energy or to create micro-electro-mechanical systems with specific properties working as sensors or actuators.

Beams are commonly used as structural elements in macro-, micro- or even nanoscales. Independently of the scale, for moderate and large vibrations nonlinear beam models have to be taken into account, and their nonlinear dynamical response must be investigated. In this respect, the frequency–response curve (FRC), which reports the response amplitude as a function of the excitation frequency, is an essential tool that summarizes the main dynamical properties of the nonlinear oscillations of the investigated system. This curve is built by harmonically forcing the vibrations of the structure, a beam in this case, up to obtaining steady-state response, and then by recording the maximum amplitude versus the frequency of the excitation.

For nonlinear systems, multiple solutions may appear at a single frequency, and their actual occurrence is related to initial conditions being in the basin of attraction of one or another solution. The FRC, which are “vertical” in the linear case (i.e. the amplitude is independent of the frequency in free oscillations), may bend towards the left (softening behaviour) or right (hardening behaviour), depending on system parameters such as the order of the natural frequency, material and geometrical nonlinearities and different boundary conditions.

The main goal of this work is to study how the FRC of the nonlinear transversal oscillations of a planar beam is affected by the boundary condition in the axial direction. With this in mind, we consider a beam which is hinged on the left end and which has a spring in the axial direction on the right end, where the transversal displacements are restrained. Varying the spring stiffness from 0 to $\infty$, we can continuously pass from the simply supported to the hinged case.

The problem of determining the nonlinear dynamics of beams, in particular their hardening/softening behaviour, has been previously addressed in the literature. For example, Atluri [1] reported the nonlinear equation of motion describing large vibrations of a hinged beam taking into account longitudinal and rotatory inertia. Using the Galerkin method, he reduced the problem to coupled ordinary differential equations and then applied the multiple timescale (MTS) method to obtain an approximate analytical solution. Silva and Glynn considered a nonplanar, flexural–torsional–inextensional beam, keeping geometrical and inertia nonlinearities up to third order [2] and, similarly to the previous author, by using a perturbation technique investigated its forced motion [3]. Later, the finite element method (FEM) has been applied to in-plane beam analysed by Reddy and Singh [4]. They tested large-amplitude free vibrations and obtained the related backbone curves of beams with various boundary conditions. The planar large-amplitude forced vibrations of clamped–clamped, clamped–hinged, hinged–hinged, hinged–simply supported and cantilever beams were evaluated by Luongo et al. [5]. They reported that beams with immovable ends (two edges are blocked in the longitudinal direction) are hardening, while systems with one end free to move in the axial direction are softening phenomenon.

Further analyses were devoted to nonlinear resonances of other systems where the beam was coupled with other elements like cables or springs. 3DOF cable-stayed beam [6], hinged–hinged beam with one end rotational spring [7], hinged–hinged and simply supported unshearable beams with one end axial spring [8] may serve as examples. Shibata et al. [9] analytically and experimentally determined a parametric resonance of hinged–simply supported beam with an axial spring with the aim of applying a passive control of FRCs by varying spring stiffnesses. Also Arumi and Yabuno [10] developed a model of a beam like system including a large mass fastened to the movable end. These authors demonstrated the important effect of the higher-order geometrical nonlinearities.

In recent years, a series of papers have been published on nonlinear beam modelling [11,12,13,14]. A shearable homogeneous beam model, with coupled axial, rotational and transverse dynamics as well as nonlinear geometrical terms resulting from large deformations, has been considered. The exact beam model with no approximations or reductions has been obtained. However, the quoted works consider only the nonlinear free oscillations. This paper is a continuation of those analyses, and here, forced nonlinear oscillations have been studied (i) by using a different analytical method and (ii) by also using numerical simulations, aimed at checking the analytical results and at highlighting other aspects of the nonlinear dynamics. Also experimental validation is planned, but it is left for future works. The presented analytical and numerical results are obtained for a very general case of thick, planar, homogeneous beams in which shear deformation and axial and rotational inertia are considered.

Fig. 1

Configuration of the initially straight beam with one end linear spring

Fig. 2

Kinematics (a) and coplanar forces (b) of the deformed beam

The paper is organized as follows. In Sect. 2, the coupled partial differential equations of motion of the Timoshenko beam are derived. The model takes into account all geometrical nonlinearities; axial, transversal and rotational inertial forces; viscous damping and transversal periodic excitation. Analytical FRCs are obtained in Sect. 3 by the MTS method, while their numerical counterparts are determined in Sect. 4 by the FEM, and they are compared in Sect. 5. Some interesting dynamical phenomena observed in the numerical simulations are reported in Sect. 6. The paper ends with conclusions and future developments (Sect. 7).

2 The beam model

We consider the Timoshenko beam shown in Fig. 1. Its configuration is described by “axial” W(Z, T) and “transversal” U(Z, T) displacements and by the rotation of the beam cross-section $\theta (Z,T)$ (Fig. 2a). Z is the spatial coordinate, which in the rest configuration denotes the distance from the left end, T, is time. The beam has initial length L, so that $0\le Z \le L$, constant cross-section A and constant second moment of inertia of the beam cross-section J. The beam is homogeneous with constant Young modulus E, shear modulus G, density $\rho$. The right-end spring stiffness is $k_s$.

The constitutive, kinematics and balance equations reported in [11] are extended by adding linear viscous damping and external excitations.

The balance of forced-damped beam element (Fig. 2b) is given by the following equations in horizontal (1), vertical (2) and rotational (3) directions, respectively:

$$\begin{aligned}&\left( N \cos { \varphi }+V \sin {\varphi }\right) ' \nonumber \\&\quad = \rho A \ddot{W}+C_W \dot{W} +P_W(Z,T), \end{aligned}$$

(1)

$$\begin{aligned}&\left( N \sin { \varphi } - V \cos {\varphi }\right) ' \nonumber \\&\quad = \rho A \ddot{U}+C_U \dot{U} + P_U(Z,T), \end{aligned}$$

(2)

$$\begin{aligned}&M'-VS'=\rho J \ddot{\theta }+C_{\theta } \dot{\theta } +P_{\theta }(Z,T), \end{aligned}$$

(3)

where the dot corresponds to the derivative with respect to time T, the prime denotes the spatial derivative with respect to axial coordinate Z, $\varphi$ is the slope angle of the beam axis, $\gamma$ is the Timoshenko angle of distortion due to shear and $S'$ is the length of deformed beam element. $C_W$, $C_U$ and $C_{\theta }$ are linear viscous damping factors in Z, X and rotational directions, respectively. M, N and V are bending moment, axial and shear forces, respectively, and their directions are shown in Fig. 2b. The external distributed loads are denoted by $P_W(Z,T)$, $P_U(Z,T)$ and $P_{\theta }(Z,T)$ (see Fig.  2b).

Assuming a linear elastic behaviour the axial, shear and bending stiffnesses are given by:

$$\begin{aligned} N=EAe, \quad V=GA\gamma , \quad M=EJk, \end{aligned}$$

(4)

where, e, the elongation of the beam axis, $\gamma$, the shear strain, and k, the curvature, are the strain measures. By simple geometrical considerations (Fig. 2a), we find:

$$\begin{aligned}&S'=\sqrt{\left( 1+W'\right) ^2+U'^2},\ \nonumber \\&\quad \cos {\varphi } =\frac{1+W'}{S'}, \ \nonumber \\&\quad \sin {\varphi } =\frac{U'}{S'}, \ \tan {\varphi }=\frac{U'}{1+W'}, \end{aligned}$$

(5)

$$\begin{aligned}&e=S'-1, \quad \gamma =\theta -\varphi , \nonumber \\&\quad k=\frac{\mathrm{d}\theta }{\mathrm{d}S}=\frac{\theta '}{S'}. \end{aligned}$$

(6)

Note that the geometrical curvature $k=\frac{\theta '}{S'}$ instead of the mechanical curvature $k=\theta '$ is used [12].

Substituting Eqs. (4)–(6) into Eqs. (1)–(3), we obtain the exact forced-damped PDEs of motion of the beam:

$$\begin{aligned}&\left\{ EA\left[ \sqrt{\left( 1+W'\right) ^2+U'^2}-1\right] \frac{1+W'}{\sqrt{\left( 1+W'\right) ^2+U'^2}} \right. \nonumber \\&\quad + \, GA \left[ \theta -\arctan {\left( \frac{U'}{1+W'}\right) }\right] \nonumber \\&\quad \left. \frac{U'}{\sqrt{\left( 1+W'\right) ^2+U'^2}}\right\} ' \nonumber \\&\quad = {\rho A \ddot{W}+C_W \dot{W}}+P_W(Z,T), \end{aligned}$$

(7)

$$\begin{aligned}&\left\{ EA\left[ \sqrt{\left( 1+W'\right) ^2+U'^2}-1\right] \right. \nonumber \\&\quad \frac{U'}{\sqrt{\left( 1+W'\right) ^2+U'^2}} + \nonumber \\&\qquad - \,GA\left[ \theta -\arctan {\left( \frac{U'}{1+W'}\right) }\right] \nonumber \\&\quad \left. \frac{1+W'}{\sqrt{\left( 1+W'\right) ^2+U'^2}}\right\} ' \nonumber \\&\quad = \rho A \ddot{U}+C_U \dot{U}+P_U(Z,T), \end{aligned}$$

(8)

$$\begin{aligned}&\quad \left[ EJ \frac{\theta '}{\sqrt{\left( 1+W'\right) ^2+U'^2}}\right] '+ \nonumber \\&\quad -\, GA\left[ \theta -\arctan {\left( \frac{U'}{1+W'}\right) }\right] \sqrt{\left( 1+W'\right) ^2+U'^2} \nonumber \\&\quad \quad = \rho J \ddot{\theta }+C_{\theta } \dot{\theta } +P_{\theta }(Z,T). \end{aligned}$$

(9)

The associated boundary conditions are:

$$\begin{aligned}&W(0,T)=0, \,\, U(0,T)=0, \,\, U(L,T)=0, \nonumber \\&\quad M(0,T)=0, \,\, M(L,T)=0, \end{aligned}$$

(10)

while the axial spring on the right-end side provides:

$$\begin{aligned}&N(L,T) \cos {\varphi } + V(L,T) \sin {\varphi } \nonumber \\&\quad + \,k_s W(L,T)=0. \end{aligned}$$

(11)

Note that for $k_s=0$ the beam becomes simply supported, while $k_s \rightarrow \infty$ implies $W(L,T)=0$ and the beam is hinged–hinged. Other cases of $k_s$ describe hinged–simply supported–spring systems. In the following analysis, the dimensionless coefficient $\kappa =\frac{k_s L}{EA}$ will be used to denote the spring stiffness.

Although the presented formulation is valid for any kind of load, in the following analysis we consider only the case of a concentrated force acting in the middle point of the beam in the X direction, harmonic in time and with amplitude $P_v$:

$$\begin{aligned}&P_W(Z,T)=0,\,\, P_U(Z,T)= P_v \delta \left( Z{-}\frac{L}{2}\right) \cos {\left( {\varOmega }T\right) }, \nonumber \\&P_{\theta }(Z,T)=0, \end{aligned}$$

(12)

where $\delta (.)$ is the Dirac delta function and ${\varOmega }$ is the frequency of excitation. Note that this load is symmetric with respect to the middle point $Z=L/2$.

Different types of load, still harmonic in time, can be easily considered by the same analytical and numerical techniques used in the following sections.

3 Analytical approach: multiple timescales method

In this section, we study analytically the nonlinear resonance phenomenon, which occurs when the frequency of excitation is near the nth linear (transversal) frequency of the beam, corresponding to a bending mode. Also sub- and superharmonic resonances (that occur when the frequency of the excitation is a integer multiple or a integer fraction of the natural frequencies) can be considered analytically in the same way, but we leave these analyses, that require a lot of computations, for further development to limit the length of the paper.

The governing PDEs of motion (7)–(9) with boundary conditions (10) and (11) are attacked directly, without introducing any reduced-order model or condensation. However, since we apply the perturbation method (PM) up to the third order, higher-order terms are not considered in this work. Thus, the PDEs of motion (7)–(9) are consistently expanded in Taylor series up to the third order:

$$\begin{aligned}&EA \left( W'+\frac{1}{2} U'^2-U'^2W '\right) ' \nonumber \\&\qquad +\,GA\left( U'\theta -U'^2+2U'^2W'-U'W'\theta \right) ' \nonumber \\&\quad ={\rho A \ddot{W}+C_W \dot{W}}, \end{aligned}$$

(13)

$$\begin{aligned}&EA\left( U'W'+\frac{1}{2}U'^3-U' W'^2\right) ' \nonumber \\&\qquad +\,GA\left( U'-\theta -U'W'+\frac{1}{2}U'^2\theta \right. \nonumber \\&\qquad \left. -\,\frac{5}{6} U'^3+U'W'^2\right) ' \nonumber \\&\quad =\rho A \ddot{U}+C_U \dot{U}+P_v \delta \left( Z-\frac{L}{2} \right) \cos {\left( {\varOmega }T\right) }, \nonumber \\ \end{aligned}$$

(14)

$$\begin{aligned}&\quad EJ\left( \theta '-W'\theta '-\frac{1}{2}U'^2\theta '\right) ' \nonumber \\&\qquad +\,GA\left( U'-\theta -W'\theta -\frac{1}{2}U'^2\theta +\frac{1}{6}U'^3\right) \nonumber \\&\quad =\rho J \ddot{\theta }+C_{\theta } \dot{\theta }. \end{aligned}$$

(15)

It must be remarked that this expansion is not an extra approximation, but it is only needed to represent in a simpler way the PM, which inherently contains the third-order approximation.

System (13)–(15) is solved analytically by the multiple timescales method [15]. The fast scale $t_0$ and two slow scales $t_1$ and $t_2$ of time are introduced:

$$\begin{aligned} t_0=T, \ t_1=\varepsilon T, \ t_2=\varepsilon ^2 T, \end{aligned}$$

(16)

where $\varepsilon$ is a formal small book-keeping parameter. Damping coefficients are assumed to be proportional to $\varepsilon ^2$, while the external excitation to $\varepsilon ^3$:

$$\begin{aligned}&C_W=\varepsilon ^2 c_W, \ C_U=\varepsilon ^2 c_U, \nonumber \\&\quad C_{\theta }=\varepsilon ^2 c_{\theta }, \ P_v=\varepsilon ^3 p_v. \end{aligned}$$

(17)

The excitation frequency is assumed to be close to a given natural frequency $\omega _n$, i.e. ${\varOmega }= \omega _n+\varepsilon ^2 \sigma$, where $\sigma$ is the detuning parameter (we should consider ${\varOmega }= \omega _n /m+\varepsilon ^2 \sigma$ for superharmonic resonance and ${\varOmega }= m \omega _n+\varepsilon ^2 \sigma$ for subharmonic resonances). Thus, we may write:

$$\begin{aligned} {\varOmega }T= & {} \left( \omega _n+\varepsilon ^2 \sigma \right) T=\left( T\omega _n+ T \varepsilon ^2 \sigma \right) \nonumber \\= & {} t_0\omega _n+t_2\sigma . \end{aligned}$$

(18)

The solutions up to third order are sought in the following forms:

$$\begin{aligned} W(Z,T)= & {} \varepsilon W_1(Z,t_0,t_1,t_2)+\varepsilon ^2 W_2(Z,t_0,t_1,t_2) \nonumber \\&+\,\varepsilon ^3 W_3(Z,t_0,t_1,t_2), \end{aligned}$$

(19)

$$\begin{aligned} U(Z,T)= & {} \varepsilon U_1(Z,t_0,t_1,t_2)+\varepsilon ^2 U_2(Z,t_0,t_1,t_2) \nonumber \\&+\,\varepsilon ^3 U_3(Z,t_0,t_1,t_2), \end{aligned}$$

(20)

$$\begin{aligned} \theta (Z,T)= & {} \varepsilon \theta _1(Z,t_0,t_1,t_2)+\varepsilon ^2 \theta _2(Z,t_0,t_1,t_2) \nonumber \\&+\,\varepsilon ^3 \theta _3(Z,t_0,t_1,t_2). \end{aligned}$$

(21)

Using the chain rule, time derivatives with respect to time T are transformed as follows:

$$\begin{aligned} \dot{W}= & {} \left( D_0 + \varepsilon D_1 + \varepsilon ^2 D_2\right) W, \end{aligned}$$

(22)

$$\begin{aligned} \ddot{W}= & {} \left[ D_0^2 + 2 \varepsilon D_0 D_1 +\varepsilon ^2 \left( 2 D_0 D_2 + D_1^2\right) \right] W,\nonumber \\ \end{aligned}$$

(23)

where $D_j$ corresponds to $\frac{\partial }{\partial t_j}$. Similar expressions are obtained for the other unknowns U and $\theta$.

Substituting expressions (19)–(23) in the governing equations and boundary conditions, and collecting terms for different powers of $\varepsilon$, we get sets of equations in successive perturbation orders.

First order

$$\begin{aligned}&EA W_1''-\rho A D_0^2W_1=0, \end{aligned}$$

(24)

$$\begin{aligned}&GA(U_1'-\theta _1)'-\rho A D_0^2 U_1=0, \end{aligned}$$

(25)

$$\begin{aligned}&EJ \theta _1''-GA (\theta _1-U_1')-\rho J D_0^2\theta _1=0. \end{aligned}$$

(26)

Boundary conditions:

$$\begin{aligned}&U_1(0,T)=0, \ U_1(L,T)=0,\ \theta _1'(0,T)=0 \nonumber \\&\theta _1'(L,T)=0,\ W_1(0,T)=0, \end{aligned}$$

(27)

$$\begin{aligned}&EA W_1'(L,T)+k_s W_1(L,T)=0. \end{aligned}$$

(28)

Second order

$$\begin{aligned}&EA W_2''-\rho A D_0^2 W_2\nonumber \\&\quad = 2 \rho A \, D_0 D_1 W_1 {-}GA \left( \theta _1 U_1' {-} U_1'^2\right) '{-}\frac{1}{2} EA \left( {U_1'^2}\right) ', \end{aligned}$$

(29)

$$\begin{aligned}&GA (U_2'-\theta _2)'-\rho A D_0^2 U_2\nonumber \\&\quad =2 \rho A D_0 D_1 U_1 -\,EA (U_1'W_1')'+GA (U_1'W_1')', \end{aligned}$$

(30)

$$\begin{aligned}&EJ \theta _2''-GA (\theta _2-U_2')-\rho J D_0^2\theta _2\nonumber \\&\quad =2\rho J D_0D_1\theta _1+\,GA (W_1' \theta _1) + EJ (W_1'\theta _1')'. \end{aligned}$$

(31)

Boundary conditions:

$$\begin{aligned}&U_2(0,T)=0, \ U_2(L,T)=0,\ \theta _2'(0,T)=0 \nonumber \\&\theta _2'(L,T)=0,\, W_2(0,T)=0, \end{aligned}$$

(32)

$$\begin{aligned}&EA W_2'(L,T)+k_s W_2(L,T)+\,\frac{1}{2} EA U_1'^2(L,T)\nonumber \\&\quad -\, GA \left[ U_1'^2(L,T) +\theta _1(L,T)U_1'(L,T)\right] =0.\nonumber \\ \end{aligned}$$

(33)

Third order

$$\begin{aligned}&EA W_3''-\rho A D_0^2 W_3=c_W D_0 W_1 \nonumber \\&\quad +\,\rho A (D_1^2+2D_0D_2)W_1 \nonumber \\&\quad +\,2\rho A D_0D_1W_2-EA (U_1'U_2'-U_1'^2W_1')'+\nonumber \\&\quad - \,GA \left( 2W_1'U_1'^2+\theta _2 U_1'-2U_2'U_1' \right. \nonumber \\&\quad \left. -\,\theta _1W_1'U_1'+\theta _1U_2'\right) ',\end{aligned}$$

(34)

$$\begin{aligned}&GA (U_3'-\theta _3)'-\rho A D_0^2 U_3 = c_U D_0 U_1 \nonumber \\&\quad +\, p_v \delta \left( Z-\frac{L}{2}\right) \cos {\left( t_0\omega _0+t_2\sigma \right) } \nonumber \\&\quad +\,\rho A (2D_0D_2+D_1^2)U_1+2\rho A D_0D_1U_2 \nonumber \\&\quad +\,EA \left( U_1'W_1'^2-U_1'W_2'-\frac{1}{2}U_1'^3-U_2'W_1' \right) ' \nonumber \\&\quad +\,GA \left( -\frac{1}{2}\theta _1U_1'^2-U_1W_1'^2+U_1'W_2' \right. \nonumber \\&\quad \left. +\,\frac{5}{6}U_1'^3+U_2'W_1' \right) ', \end{aligned}$$

(35)

$$\begin{aligned}&EJ \theta _3''-GA\left( \theta _3-U_3'\right) -\rho J D_0^2\theta _3 \nonumber \\&\quad =c_{\theta } D_0\theta _1+\rho J \left( 2D_0D_2+D_1^2\right) \theta _1 \nonumber \\&\qquad +\,\rho J 2D_0D_1\theta _2 \nonumber \\&\qquad +\,GA\left( \frac{1}{2}\theta _1U_1'^2-\frac{1}{6}U_1'^3+\theta _2W_1'+\theta _1W_2'\right) \nonumber \\&\qquad +\,EJ\left( \frac{1}{2}\theta _1'U_1'^2-\theta _1'W_1'^2+\theta _2'W_1'+\theta _1'W_2'\right) '. \end{aligned}$$

(36)

Boundary conditions:

$$\begin{aligned}&U_3(0,T)=0, \ U_3(L,T)=0,\ \theta _3'(0,T)=0 \nonumber \\&\quad \theta _3'(L,T)=0,\ W_3(0,T)=0, \end{aligned}$$

(37)

$$\begin{aligned}&EA W_3'(L,T)+k_s W_3(L,T) \nonumber \\&\quad +\,GA [-\theta _1(L,T)W_1'(L,T)U_1'(L,T) \nonumber \\&\quad +\,\theta _2(L,T)U_1'(L,T)+U_2'(L,T)\theta _1(L,T) ] \nonumber \\&\quad +\,(EA-2GA) \nonumber \\&\quad \left[ U_2'(L,T)U_1'(L,T) \right. \nonumber \\&\quad \left. -\,W_1'(L,T)U_1'^2(L,T) \right] =0. \end{aligned}$$

(38)

3.1 First-order solution

The first-order solution is constituted by two uncoupled eigenvalue problems, one in the axial direction, involving $W_1$, and the other in the transversal direction, involving $U_1$ and $\theta _1$. As in [12], and differently from [13], in this paper we consider only transversal behaviour to the first order; thus,

$$\begin{aligned} W_1(Z,T)=0. \end{aligned}$$

(39)

The remaining part of the solution is given by:

$$\begin{aligned}&U_1(Z,t_0,t_1,t_2)=\left[ A_{\mathrm{re}}(t_1,t_2) \mathrm{e}^{i \omega _n t_0}+A_{{rm im}}(t_1,t_2) \right. \nonumber \\&\quad \left. \mathrm{e}^{-i \omega _n t_0}\right] \hat{U}_{1,n}(Z), \end{aligned}$$

(40)

$$\begin{aligned}&\theta _1(Z,t_0,t_1,t_2) \nonumber \\&\quad =\left[ A_{\mathrm{re}}(t_1,t_2) \mathrm{e}^{i \omega _n t_0} \right. \nonumber \\&\qquad \left. +\,A_{\mathrm{im}}(t_1,t_2)\mathrm{e}^{-i \omega _n t_0}\right] \hat{\theta }_{1,n}(Z), \end{aligned}$$

(41)

where $\hat{U}_{1,n}(Z)$ and $\hat{\theta }_{1,n}(Z)$ are the nth modal shapes:

$$\begin{aligned}&\hat{U}_{1,n}(Z)=\sin {\left( \frac{n\pi Z}{L}\right) }, \nonumber \\&\hat{\theta }_{1,n}(Z)=\left( \frac{n \pi }{L}-\frac{L \rho \omega _n^2}{\pi G n} \right) \cos {\left( \frac{n \pi Z}{L}\right) }. \end{aligned}$$

(42)

The nth natural frequency of the beam is given by

$$\begin{aligned}&\omega _n^2=\frac{A G}{2 J \rho }+\frac{\pi ^2 n^2 (E+G)}{2 L^2 \rho } \nonumber \\&\quad -\,\frac{\sqrt{\left[ A G L^2 +\pi ^2 J n^2(E+G) \right] ^2-4 \pi ^4 E G J^2 n^4}}{2 J L^2 \rho }.\nonumber \\ \end{aligned}$$

(43)

The above solution means that we are considering only a single primary resonances. If we want to investigate also internal resonances, in (40) and (41) we have to consider two terms with two natural frequencies $\omega _{n}$ and $\omega _{m}$ where $\omega _{n}/\omega _{m}$ or $\omega _{m}/\omega _{n}$ is close to a natural number.

3.2 Second-order solution

The solution of second-order Eq. (29) has the form

$$\begin{aligned} W_2(Z,t_0,t_1,t_2)= & {} W_{2a}(Z,t_1,t_2) \nonumber \\&\quad +\,W_{2b}(Z,t_1,t_2)\mathrm{e}^{2i \omega _n t_0} \nonumber \\&\quad +\,W_{2c}(Z,t_1,t_2) \mathrm{e}^{-2i \omega _n t_0}, \end{aligned}$$

(44)

where

$$\begin{aligned} W_{2a}(Z,t_1,t_2)= & {} -A_{\mathrm{im}} \left( t_1,t_2\right) A_{\mathrm{re}} \left( t_1,t_2\right) \left( E n^2 \pi ^2-2 L^2 \rho \omega _n^2 \right) \nonumber \\&\quad \times \, \frac{2 n \pi Z+L\left( \kappa +1\right) \sin {\left( \frac{2 n \pi Z}{L}\right) }}{4 E L^2 n \pi \left( \kappa +1\right) }, \end{aligned}$$

(45)

$$\begin{aligned}&W_{2b}(Z,t_1,t_2) \nonumber \\&= - \frac{A_{\mathrm{re}}\left( t_ 1,t_ 2\right) ^2 \left( \pi ^2 E n^2 -2 L^2 \rho \omega _n^2\right) }{\sqrt{E}\kappa \sin {\left( \frac{2 L \sqrt{\rho } \omega _n }{\sqrt{E}}\right) }+2 L \sqrt{\rho } \omega _n \cos {\left( \frac{2 L \sqrt{\rho }\omega _n}{\sqrt{E}}\right) }} \nonumber \\&\qquad \times \, \frac{1}{16 \sqrt{E} L \left( E n^2 \pi ^2-L^2 \rho \omega _n^2\right) } \nonumber \\&\qquad \times \,\left[ 4 \pi ^2 E n^2 \sin \left( \frac{2 \sqrt{\rho } Z \omega _n}{\sqrt{E}}\right) \right. \nonumber \\&\qquad +\, 2\pi E \kappa n \sin \left( \frac{2 \pi n Z}{L}\right) \sin \left( \frac{2 L \sqrt{\rho } \omega _n}{\sqrt{E}}\right) \nonumber \\&\qquad -\,8 \rho \omega _n^2 L^2 \sin \left( \frac{2 \sqrt{\rho } Z \omega _n}{\sqrt{E}}\right) \nonumber \\&\qquad +\,4 \pi \sqrt{E} L n \sqrt{\rho } \omega _n \sin \left( \frac{2 \pi n Z}{L}\right) \nonumber \\&\qquad \left. \cos \left( \frac{2 L \sqrt{\rho } \omega _n}{\sqrt{E}}\right) \right] , \end{aligned}$$

(46)

$$\begin{aligned}&W_{2c}(Z,t_1,t_2) \nonumber \\&\quad =\frac{A_{\mathrm{im}}\left( t_ 1,t_ 2\right) ^2}{A_{\mathrm{re}}\left( t_ 1,t_ 2\right) ^2}W_{2b}(Z,t_1,t_2). \end{aligned}$$

(47)

Second-order Eqs. (30) and (31) and the associated boundary conditions (32) and (33), considering assumption (39), can be solved if and only if the following solvability condition is satisfied:

$$\begin{aligned}&\int _{0}^{L} \left( 2\rho A D_0D_1U_1\right) U_1 \mathrm{d}Z \nonumber \\&\quad +\, \int _{0}^{L} \left( 2\rho J D_0D_1\theta _1\right) \theta _1 \mathrm{d}Z=0. \end{aligned}$$

(48)

It provides

$$\begin{aligned} \frac{\partial A_{\mathrm{re}}}{\partial t_1}=0, \ \frac{\partial A_{\mathrm{im}}}{\partial t_1}=0; \end{aligned}$$

(49)

namely, the real and imaginary amplitudes are independent of the slow time $t_1$, $A_{\mathrm{re}}(t_2)$ and $A_{\mathrm{im}}(t_2)$.

With (49) and (39), system (30) and (31) simplifies to

$$\begin{aligned}&GA (\theta _2-U_2')'+\rho A D_0^2 U_2=0, \end{aligned}$$

(50)

$$\begin{aligned}&EJ \theta _2''-GA (\theta _2-U_2')-\rho J D_0^2\theta _2=0, \end{aligned}$$

(51)

and its solution is given by

$$\begin{aligned}&U_2(Z,t_0,t_2)= 0, \end{aligned}$$

(52)

$$\begin{aligned}&\theta _2(Z,t_0,t_2)=0. \end{aligned}$$

(53)

3.3 Third-order solution

In the third-order approximation, we do not need the approximate solution, but just the solvability condition

$$\begin{aligned}&\int _{0}^{L} \left[ \rho A 2D_0D_2 U_1 \right. \nonumber \\&\quad + \,GA \left( -\frac{1}{2}\theta _1U_1'^2+U_1'W_2'+\frac{5}{6}U_1'^3 \right) ' \nonumber \\&\quad +\, EA \left( -U_1'W_2'-\frac{1}{2}U_1'^3\right) '+ c_U D_0 U_1 \nonumber \\&\quad +\, \left. p_v \delta \left( Z- \frac{L}{2} \right) \cos {\left( t_0\omega _n+t_2\sigma \right) }\right] U_1 \mathrm{d}Z \nonumber \\&\quad +\,\int _{0}^{L} \left[ \rho J 2D_0D_2\theta _1+GA \right. \nonumber \\&\quad \left. \left( \frac{1}{2}\theta _1U_1'^2-\frac{1}{6}U_1'^3+\theta _1W_2'\right) \right. \nonumber \\&\quad +\,\left. EJ\left( \frac{1}{2}\theta _1'U_1'^2+\theta _1'W_2'\right) ' +c_{\theta }D_0\theta _1\right] \theta _1 \mathrm{d}Z =0\nonumber \\ \end{aligned}$$

(54)

of Eqs. (35) and (36) and boundary conditions (37) and (38) is enough for our purposes.

After cumbersome computations, separating the terms containing $\mathrm{e}^{i\omega _n t_0}$ from those containing $\mathrm{e}^{-i\omega _n t_0}$, Eq. (54) gives the following two ordinary differential equations in the unknowns $A_{\mathrm{re}}$ and $A_{\mathrm{im}}$:

$$\begin{aligned}&i c_1 \frac{\partial A_{\mathrm{re}}}{\partial t_2} + i c_2 A_{\mathrm{re}} + 4 c_3 A_{\mathrm{re}}^2A_{\mathrm{im}} \nonumber \\&\quad + \,\frac{1}{2} p_v \mathrm{e}^{ i \sigma t_2} \sin \left( \frac{n \pi }{2}\right) = 0,\end{aligned}$$

(55)

$$\begin{aligned}&-\, i c_1 \frac{\partial A_{\mathrm{im}}}{\partial t_2} - i c_2 A_{\mathrm{im}} + 4 c_3 A_{\mathrm{im}}^2A_{\mathrm{re}} \nonumber \\&\quad +\, \frac{1}{2} p_v \mathrm{e}^{ -i \sigma t_2} \sin \left( \frac{n \pi }{2}\right) = 0, \end{aligned}$$

(56)

where

$$\begin{aligned} c_1= & {} \rho \omega _n \frac{J\left( G n^2 \pi ^2 - L^2 \rho \omega _n^2\right) ^2+AL}{ G^2 L n^2 \pi ^2 }, \nonumber \\ c_2= & {} \omega _n \left[ c_\theta \frac{ \left( G n^2 \pi ^2 -L^2 \rho \omega _n^2\right) ^2}{2 G^2 L n^2 \pi ^2}+c_U \frac{L}{2} \right] , \end{aligned}$$

(57)

and where the real value coefficient $c_3$, has a very long expression that cannot be reported directly but can be handled by a symbolic manipulator software.

The coefficient $c_1$ takes into account inertia, while $c_2$ damping coefficients, and they both have positive values. The key coefficient $c_3$, on the other hand, summarizes all the nonlinear effects and depends on A, J, E, G, L, n, $\kappa$, $\rho$. For illustrative purposes, it is reported in Fig. 3, where we show its dependence on $\kappa$ since, as said in Introduction, in the present work we are mainly interested in determining the effect of the spring stiffness on the nonlinear dynamical behaviour.

Fig. 3

The coefficient $c_3(\kappa )$ for $A=25~\hbox {cm}^2$, $J=52.08~\hbox {cm}^4$, $E=21{,}000~\hbox {kN}/\hbox {cm}^2$, $L=50~\hbox {cm}$, $G=8000~\hbox {kN}/\hbox {cm}^2$, $\rho =0.00785~\hbox {kg}/\hbox {cm}^3$ and for modes $n=1$, $n=2$, $n=3$

Fig. 4

The coefficient $c_3(\kappa )$ for $A=25~\hbox {cm}^2$, $J=52.08~\hbox {cm}^4$, $E=21{,}000~\hbox {kN}/\hbox {cm}^2$, $L=100~\hbox {cm}$, $G=8000~\hbox {kN}/\hbox {cm}^2$, $\rho =0.00785~\hbox {kg}/\hbox {cm}^3$ and for modes $n=1$, $n=2$, $n=3$

For increasing spring stiffness, the $c_3$ coefficient monotonically (for $n=1$ and $n=3$) tends to the value corresponding to the hinged–hinged beam. It quickly approaches the asymptote, and at $\kappa \cong 15$, the difference with the asymptotic limit is negligible.

At about $\kappa \cong 1.43$, and for the second mode ($n=2$), $c_3$ has a singularity and tends to $\pm \infty$, as a consequence of a denominator that vanishes. This singularity will be also shown in forthcoming Fig. 5c.

It is worth to remark that for even values of n the excitation terms disappear in (55) and (56) according to the fact that the given load is symmetric, while the even linear modes are antisymmetric with respect to the beam middle point.

In order to illustrate also the effect of the slenderness of the beam, in Fig. 4, we draw the function $c_3(\kappa )$ for a different length, while the other parameters are the same of Fig. 3.

In the considered range of $\kappa$, the beam length caused the shift of the singularity from the second to third mode, and the critical spring stiffness is now about $\kappa =3.21$. The asymptotic values for $\kappa \rightarrow \infty$ are slightly smaller than those for $L=50$ cm, showing how the slenderness tends to reduce the nonlinear effects for the hinged–hinged beam. The coefficient $c_3$ gets negative values for $\kappa <0.04$ ($n=1$) and $\kappa <0.2$ ($n=3$), while for higher spring stiffnesses it is positive.

3.4 Analytical frequency–response curves

To determine the frequency–response curves, we have to solve Eqs. (55) and (56). In this respect, it is useful to express the complex amplitude in the polar form

$$\begin{aligned} A_{\mathrm{re}}(t_2)=\frac{1}{2} a(t_2) \mathrm{e}^{i \beta (t_2)}, \ A_{\mathrm{im}}(t_2)=\frac{1}{2} a(t_2) \mathrm{e}^{-i \beta (t_2)}. \end{aligned}$$

(58)

From (55) and (56), we obtain

$$\begin{aligned} \frac{\partial a}{\partial t_2}= - \frac{c_2}{c_1} a - \frac{p_v}{ c_1} \sin \left( \frac{n \pi }{2}\right) \sin {(\sigma t_2 -\beta )}, \end{aligned}$$

(59)

and

$$\begin{aligned} a \frac{\partial \beta }{\partial t_2}=\frac{ c_3 }{ c_1 } a^3 + \frac{p_v}{c_1} \sin \left( \frac{n \pi }{2}\right) \cos {(\sigma t_2 -\beta )}. \end{aligned}$$

(60)

Introducing the new variable

$$\begin{aligned} \gamma (t_2)=\sigma t_2-\beta (t_2) \quad \Rightarrow \quad \frac{\partial \gamma (t)}{\partial t_2}=\sigma -\frac{\partial \beta }{\partial t_2}, \end{aligned}$$

(61)

system (59) and (60) becomes

$$\begin{aligned}&\frac{\partial a}{\partial t_2}= - \frac{c_2}{c_1} a - \frac{p_v}{ c_1} \sin \left( \frac{n \pi }{2}\right) \sin {(\gamma )}, \end{aligned}$$

(62)

$$\begin{aligned}&a \frac{\partial \gamma }{\partial t_2}=\sigma a + \frac{ c_3 }{c_1 } a^3 + \frac{p_v}{c_1} \sin \left( \frac{n \pi }{2}\right) \cos {(\gamma )}. \end{aligned}$$

(63)

Looking for the steady-state solution, we assume $\frac{\partial \gamma (t)}{\partial t_2}=0$ and $\frac{\partial a(t)}{\partial t_2}=0$. From Eqs. (62) and (63), we then obtain

$$\begin{aligned}&\sin {(\gamma )} = - \frac{c_2 a }{p_v \sin \left( \frac{n \pi }{2}\right) }, \nonumber \\&\cos {(\gamma )} = - \frac{ a \left( c_3 a^2 - c_1 \sigma \right) }{p_v \sin \left( \frac{n \pi }{2}\right) }. \end{aligned}$$

(64)

Considering that $\sin ^2{\left( \gamma \right) }+\cos ^2{\left( \gamma \right) }=1$, we get the frequency–response equation (FRE) of the nth bending mode:

$$\begin{aligned} p_v^2\sin ^2\left( \frac{n \pi }{2}\right) =c_2^2 a^2 + a^2 \left( c_3 a^2 - c_1 \sigma \right) ^2. \end{aligned}$$

(65)

Fig. 5

Function $c(\kappa )$ for first (a, b), second (c) and third (d) modes. Beam parameters as in Fig. 3

The solution of Eq. (65) provides the searched FRC, which actually is the function $a(\sigma )$. Only the inverse function $\sigma (a)$ can be expressed in a relatively simple way:

$$\begin{aligned} \sigma = \frac{c_3 a^2 \pm \sqrt{\frac{p_v^2}{a^2} \sin ^2\left( \frac{n \pi }{2}\right) - c_2^2}}{c_1} . \end{aligned}$$

(66)

On the basis of the analytical solutions, some examples will be drawn and the results will be compared with their numerical counterparts.

In the absence of load ($p_v=0$) and damping ($c_U=0$, $c_\theta =0$, i.e. $c_2=0$), expression (66) becomes

$$\begin{aligned} \sigma =c a^2, \quad c=\frac{c_3}{c_1}, \end{aligned}$$

(67)

which is known as the backbone curve and summarizes the main feature of the free nonlinear oscillations of the considered system. In particular, if the nonlinear correction coefficient c is positive, we have the hardening behaviour, and the nonlinear natural frequency increases with the amplitude of the oscillation a. On the contrary, when c is negative, the system demonstrates the softening phenomenon, and the nonlinear frequency decreases while increasing a. In the singular case $c=0$, the beam behaves as linear (up to the third-order approximation).

Table 1 First four linear natural frequencies and corresponding free vibration mode shapes

Assuming the dimensions of the beam used in Fig.  3, the coefficient c as a function of the spring stiffness $\kappa$ is depicted in Fig. 5.

We note that the first mode is softening up to $\kappa \cong 0.18$ and hardening above this threshold (see Fig. 5b). Thus, in particular, we note that the hinged–simply supported beam ($\kappa =0$) has softening behaviour, while the hinged–hinged beam ($\kappa \rightarrow \infty$) is hardening, confirming the results of the literature (see, for example, [12]).

The second mode is mainly the hardening ($c>0$), apart from a small interval around $\kappa \simeq 1.8$, where further analyses are required to better understand this particular behaviour. However, this is out of the scope of the present work and is left for further developments.

The third mode has always hardening nature. It is worth to note that the asymptotic value of c for $\kappa \rightarrow \infty$ increases by increasing n, which means that the backbone curve of higher- order modes is more and more bent towards high frequencies.

It is interesting to underline how, for a fixed value of $\kappa$, different modes can have different mechanical behaviour. For example, for $\kappa <0.18$ the first mode is softening and the second and third modes are hardening.

4 Numerical approach: Finite element method

A finite element model is created for a Timoshenko beam with initial length $L=500$ mm, cross-section $A=50 \times 50$ mm (length to width ratio equal 10), Young modulus $E=210$ GPa, Poisson’s ratio $\nu =0.3$, density $\rho =7850$ kg/$\hbox {m}^3$, beam shear factor $\chi =0.85$ and viscous damping factor $\zeta =0.06$. Use is made of the commercial software Abaqus_CAE©.

Fig. 6

First linear bending mode shape: a displacement and b rotation. Second linear bending mode shape: c displacement and d rotation. Third linear bending mode shape: e displacement and f rotation

The B31-type beam elements, having three translational degrees of freedom at each node, are used. The beam’s length is divided into 100 equal elements, which have linear shape between nodes, and allows transverse shear deformation [16]. Due to large amplitudes of vibrations, the “NLgeom” framework is applied to calculate true stress and strain instead of nominal stress and strain. It is worth to highlight that it gives valuable numerical results despite that the simplest beam element is used. Out-of-plane displacement of the nodes has been constrained to guarantee planar motion. Boundary conditions are defined as in (10) and (11). An axial linear spring of stiffness $k_s$ is inserted to link the movable end and the ground. In the case of modelling the hinged–simply supported beam, the spring has been removed, while for the hinged–hinged beam the displacement of the right beam end has been axially restrained, $W(L,T)=0$.

Fig. 7

Transversal (a, c, e) and longitudinal (b, d, f) vibrations of the hinged–simply supported-spring beam ($\kappa =4$). Global sweep $\frac{{\varOmega }}{\omega _{1}} \$increased $0\rightarrow 2$, reset and decreased $2\rightarrow 0 \$ (a, b), attractors for $\frac{{\varOmega }}{\omega _{1}}$ 1.25, 1.275, 1.3, 1.325, 1.35 (c, d), and the steady-state solution for $\frac{{\varOmega }}{\omega _{1}} \ =1.125$ (e, f)

Initially, the numerical modal analysis of the structure has been performed using Lanczos method. The natural frequencies and corresponding mode shapes of free vibrations of the beam with various spring stiffnesses have been determined. As deduced from the analytical developments, for small amplitudes, linear bending modes in the transversal direction are independent of axial spring stiffness. Furthermore, longitudinal vibrations have no effects on linear bending modes and natural frequencies.

The outcome of the linear modal analysis is reported in Table  1 and in Fig. 6. For all inspected stiffnesses $\left( 0\le k_s \le \infty \right)$, the first linear frequency is equal $\omega _1=461.47$ Hz (2899.50 $\frac{\mathrm{rad}}{\mathrm{s}}$) with period $T=0.002167$ s and corresponds to a bending mode. These numerical results perfectly match their analytical counterparts obtained in Sect. 3.1.

To study nonlinear forced vibrations, the numerical computations are performed in a dynamic explicit module. A special technique, similar to the continuation method, is adopted with the excitation frequency varied gradually. At each excitation frequency, transient computations are run for a required number of periods to reach the steady-state solution as presented in Fig. 7c, d. The achieved amplitude is recorded (Fig. 7e, f), and then, the frequency slightly changed. New numerical simulations start from the initial conditions of the steady state of the previous case, so that it is expected that the transient is short, reducing computational time and allowing to follow a given path of the solution without an unwanted jump due to the multistability and bad choice of the initial conditions. This procedure enables building the frequency–response curve “easily” and almost automatically—although the whole computations needed for obtaining a single curve are still time-consuming.

The amplitude of the concentrated time periodic force applied in the beam midpoint [see Eq. (12)] is $P_v=40799.2$ N, corresponding to 1 mm static deflection ($2\%$ of the beam thickness) in $Z=\frac{L}{2}$ of the hinged–simply supported beam.

The complete frequency–response curve is constructed in the range from 0 Hz up to the double of the natural frequency ($2 \times 461.47=922.94$ Hz for the first mode) by sweeping excitation frequency forward and backward, starting from static deflection (see Fig. 4a, b). The time step has been fixed as constant for each frequency and equal to 1 / 40 of the excitation period.

Fig. 8

MTS method (solid lines) versus Lindstedt–Poincaré (LP) method (dashed) backbone curves for beam material characteristics: $L=5$ m, $H=0.5$ m, $B=0.2$ m, $\rho = 7850\frac{kg}{m^3}$, $\nu =0.3$, ${E=2.11 \ 10^{11}\frac{N}{m^2}}$, shear correction factor $\chi =1.2$ [14]

Examples of numerical frequency–response curves for the considered case, similar to those obtained in [17], are reported in forthcoming Fig. 9, where they are compared with their analytical counterparts.

It is worth to remark that only the principal resonance zone with resonant and non-resonant paths has been considered in this work. Other solutions, or detached resonance curves [18,19,20], can be obtained by varying initial conditions (initial mode shape) for each frequency of excitation or by using hints from analytical investigations or experiments.

5 Comparisons

5.1 Free nonlinear oscillations

In this section, we compare results obtained by the multiple time scales method with those obtained in [12, 14] by the Lindstedt–Poincaré method. The comparison is made for planar Timoshenko beams with no damping and no forcing, and thus, only backbones curves for various spring stiffnesses are computed. The comparison between the two different perturbation methods shows an excellent agreement (Fig. 8), and the differences are negligible and occur only for very large values of the amplitude.

5.2 Forced nonlinear oscillations

The frequency–response curves obtained in Sect. 5 are illustrated in Fig. 9, where they are compared with their analytical counterparts obtained in Sect. 4.

The hinged–simply supported beam demonstrates softening FRC, as the upper branches of the curve bends towards the left, i.e. the frequency of oscillation decreases, while amplitude of beam response increases.

Fig. 9

The primary resonance frequency–response curves of hinged–simply supported–spring systems obtained by FEM (dots) and MTS method (continuous lines) for selected spring stiffnesses. Characteristic changes from softening to hardening, from the left to right, respectively

By increasing the stiffness of the end spring, the qualitative behaviour changes. In particular, we note that even for relatively small values of $\kappa$ the beam becomes hardening, and the curve bends towards the right. This effect increases while increasing $\kappa$, up to the limit case $\kappa =\infty$ corresponding to the hinged–hinged boundary conditions.

As it can be observed in Fig. 9, a good agreement is obtained between analytical and numerical simulations. It is worth to remark that, for small values of $\kappa$, this agreement is also quantitatively excellent, even for large amplitudes. For moderate and large values of $\kappa$, we still have a very good quantitative agreement up to moderate amplitudes, about 0.02 m. Above this threshold, there are some discrepancies between numerical and analytical results. For example, for $\sigma =\omega _n \ (2899.5\,\frac{\mathrm{rad}}{\mathrm{s}})$ the amplitudes obtained by FEM are $27.74\%$ ($\kappa =\infty$), $26.87\%$ ($\kappa =4$) and $28.55\%$ ($\kappa =1$) higher than those obtained analytically by the MTS method. This difference is likely due to the fact that for large amplitude the MTS method is not expected to be very precise, as in this paper it is developed up to the third order. Going up to higher perturbation order likely, the accuracy of the analytical solution will be improved. Actually, the discrepancies for large values of $\kappa$ and large values of the amplitudes are not surprising, but rather it is positively surprising there is a numerical/analytical agreement that for low values of $\kappa$ and very large amplitudes (0.07 m, i.e. $14\%$ of the length of the beam), much above than expected. The above results represent a cross-check of the reliability of both numerical and analytical results.

Fig. 10

Fast Fourier transform of steady-state vibrations of beam with spring stiffness $\kappa =4$ for $\frac{{\varOmega }}{\omega _{1}}=1.4$. Transformation of time history at $Z=0.75L$ in the U-direction

Fig. 11

a Time histories of middle point of the beam in longitudinal direction for $\kappa =1$, excitation frequencies $\frac{{\varOmega }}{\omega _{1}}=1.2$, $\frac{{\varOmega }}{\omega _{1}}=1.225$, $\frac{{\varOmega }}{\omega _{1}}=1.25$ from the left to right, respectively. b Double period of steady oscillations $\frac{{\varOmega }}{\omega _{1}}=1.225$. Fast Fourier transform of steady-state vibrations for c $\frac{{\varOmega }}{\omega _{1}}=1.2$ and d $\frac{{\varOmega }}{\omega _{1}}=1.225$

For hinged–simply supported beams with axial spring stiffnesses $\kappa =\frac{1}{4}$ and $\kappa =\frac{1}{8}$, premature jumps from upper to lower brunches occur for $\frac{{\varOmega }}{\omega _{1}}=1.275$ and $\frac{{\varOmega }}{\omega _{1}}=1.175$, respectively. Despite small frequency growth the initial conditions (upshot from previous frequency steady-state solution), go outside the basin of attraction and achieve the solution in the lower path. This phenomenon can be a result of very small and eroded basins of attraction, which even for the minor perturbation of frequency produces escape from the basins. A detailed analysis of this interesting aspect requires the use of dynamical integrity ideas and tools [21, 22], but this is out of the scope of the present work.

6 Superharmonic and internal resonances

As previously said, analytical approach is limited to individual nth primary resonances, without considering superharmonic/subharmonic resonances and with no interaction between nth and mth modes or longitudinal and transversal harmonic components. However, these phenomena appear naturally in numerical simulations, and thus, in this section numerical results based on FEM simulations are inspected to preliminary check the existence of possible superharmonic resonances and the multimode engagement. The detailed analytical study of these phenomena is left for future works.

Superharmonic resonances appear for spring stiffnesses $\kappa =1$ and $\kappa =4$ (see Fig. 9), in the neighbourhood of the frequencies $\frac{{\varOmega }}{\omega _{1}}=1.75$ and $\frac{{\varOmega }}{\omega _{1}}=1.4$, respectively. They increase the amplitude of about $8.897\%$ ($\kappa =4$) and $2.84\%$ ($\kappa =1$) with respect to the trend of backbone curves. This phenomenon may generate a new solution path. A deep investigation of existence of possible new solutions is out of the scopes of the present paper and is left for future work.

The fast Fourier transform (FFT) of Fig. 10 demonstrates that for $\kappa =4$, in addition to the main peak at $\frac{\omega }{\omega _1}=1.4$, i.e. $\omega ={\varOmega }$, there is a second peak at $\frac{\omega }{\omega _1}=4.2$, i.e. $\omega =3 {\varOmega }$. The FFT of the other case, not reported here to limit the length of the paper, shows that for $\kappa =1$ there is a main peak at $\frac{\omega }{\omega _1}=1.725$ and a second peak at $\frac{\omega }{\omega _1}=5.175$. The ratio in both cases is equal 3 showing the occurrence of a superharmonic resonance. Those peaks correspond to the first and second bending modes, despite that the excitation is applied to middle point of the beam in transversal direction, so that it directly excites the first bending mode, while it does not trigger (directly) the second bending mode, which is antisymmetric and has a nodal point for $Z=\frac{L}{2}$. Thus, this effect clearly occurs due to the nonlinear coupling between the first and the second modes and represents a nonlinear internal resonance.

Another interesting phenomenon is an axial–transversal internal resonance obtained for spring stiffness $\kappa =1$ and corresponding to frequencies $1.2\le \frac{{\varOmega }}{\omega _{1}}\le 1.25$ of the increasing sweep. In this case, to highlight the involved dynamical phenomenon, we used a time step equal $\frac{1}{160}$ of the excitation period instead of the current $\frac{1}{40}$ (Fig.  11a, b). To better understand the result of Fig. 11, we remind that, due to the symmetry of the system with respect to the rest position, the period of longitudinal vibrations is half of the period of transverse vibrations, and the frequency is double.

The FFT of the axial vibrations of a centre point of the beam shows that vibrations are composed of three harmonics $\omega _{11}=2{\varOmega }_{1.2}$, $\omega _{12}=4{\varOmega }_{1.2}$, $\omega _{13}=6{\varOmega }_{1.2}$ (Fig. 11c) and $\omega _{21}=2{\varOmega }_{1.225}$, $\omega _{22}=4{\varOmega }_{1.225}$, $\omega _{23}=6{\varOmega }_{1.225}$ (Fig. 11d). Frequencies $\omega _{11}$, $\omega _{12}$, $\omega _{13}$ and $\omega _{21}$, $\omega _{21}$ correspond to longitudinal vibrations accompanying the first bending mode. In this case, the internal resonance between axial and transversal vibrations is highlighted by the “large” amplitude corresponding to $\omega _{23}=3394.5\,\mathrm{Hz}$, which happens because the first longitudinal mode is excited, since its natural frequency 3340 Hz is very close to $\omega _{23}$.

7 Conclusions and further developments

The complete model of a shearable homogeneous beam taking into account coupled axial, rotational and transverse vibrations has been presented in the paper. The nonlinear geometrical terms resulting from large deformations and a nonlinear curvature have been considered in the model. The exact equations of motion of a planar, initially straight, beam including forcing and damping terms have been derived and then expanded up to third order. The hardening/softening behaviour of the hinged–simply supported–spring beam and the influence of end spring stiffnesses on frequency response curves have been determined by two methods: analytical multiple timescales and explicit numerical by finite elements.

The resonances analysis has been performed by applying the MTS method to the nth bending mode. Using the solvability conditions, the modulation equations have been obtained, and then, the frequency response curves and backbone curves have been determined.

Numerical analysis has disclosed the level of complexity of frequency–response curves and confirmed the correctness of the results of the analytical computations. It has been detected that for some cases, time histories of steady-state solutions are strongly composed of more than one harmonic. Superharmonic resonance involving the coupling between first and second bending modes significantly affected the main path of transverse vibrations. The internal resonance between transversal and longitudinal modes has been reported, too.

Many further developments are possible and worthy of future investigations. Among them, we quote:

• performing experimental tests;

• analytical solutions of subharmonic, superharmonic and internal resonances;

• study of different configurations, including different boundary conditions;

• investigation in detail the singular behaviour of the nonlinear correction coefficient domains, where it tends to infinity;

• study of the softening versus hardening dichotomy which is observed by increasing stiffness of the axial spring;

• exploitation of axial/transversal couplings;

• considering multifrequency excitations;

• investigating a nonplanar beam with the twisting effect.