Longitudinal–transversal internal resonances in Timoshenko beams with an axial elastic boundary condition

The internal resonances between the longitudinal and transversal oscillations of a forced Timoshenko beam with an axial end spring are studied in depth. In the linear regime, the loci of occurrence of 1 : ir, ir∈N\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ir \in \mathbb {N}$$\end{document}, internal resonances in the parameters space are identified. Then, by means of the multiple time scales method, the 1 : 2 case is investigated in the nonlinear regime, and the frequency response functions and backbone curves are obtained analytically, and investigated thoroughly. They are also compared with finite element numerical simulations, to prove their reliability. Attention is paid to the system response obtained by varying the stiffness of the end spring, and it is shown that the nonlinear behaviour instantaneously jumps from hardening to softening by crossing the exact internal resonance value, in contrast to the singular (i.e. tending to infinity) behaviour of the nonlinear correction coefficient previously observed (without properly taking the internal resonance into account).


Introduction
The dependence of the nonlinear transversal behaviour of beams on different axial boundary conditions is an interesting topic that has been studied by various authors. In [1], it was shown that a beam is hardening if the axial displacement is constrained at the boundary and is softening if the boundary is free to move axially. These findings were confirmed in [2], where a refined reduced order model analysis was performed, and the effect of the axial inertia and inextensibility was investigated.
A fully continuous model was considered in [3], where the multiple time scales method is applied directly to the partial differential equations of motion. The results were also compared with some experi-ments. On the same line of investigation, [4] considered an axial (i.e. parametric) excitation, and developed the asymptotic analysis up to the fifth order, i.e. one order more than what is usually done. Also some experimental results were reported.
In fact, it is known that axial inertia contributes to softening, while the geometrical stiffness, which is strongly affected by the axial constraint, contributes to hardening. The whole behaviour depends on which one of these two opposite phenomena predominates. Thus, it is possible to foresee that introducing at the boundary a spring, of stiffness κ, which allows us to regulate the axial constraint and thus the geometrical nonlinearity, permits to control the hardening/softening behaviour, passing from softening when κ = 0 (no constraint) to hardening when κ → ∞ (full constraint), that are the limit cases reported in the earlier literature [1,2].
Although a boundary axial spring is present in [3] and [4], its effect on the softening/hardening behaviour is not fully investigated there. This has been done, instead, in [5], where the dependence of the nonlinear correction coefficient ω 2 (ω 2 > 0 gives hardening, ω 2 < 0 provides softening) upon κ is deeply investigated. There, also the effects of the axial and rotational inertia have been studied, as well as those of the shear stiffness and slenderness. A Timoshenko beam model was considered to have reliable results also for thick beams.
Starting from [6], where the equations of motions were firstly derived and initially analyzed by means of the Poincaré-Lindstedt method, and [5], the authors have done many studies on this topic. In [7], the approximate analytical results have been compared with numerical simulations to check their reliability. In these papers, attention was focused on the transversal oscillations, and the axial vibrations are set equal to zero to the first order and appear only to the second order because of the nonlinear coupling. Actually, this hypothesis has been relaxed in [8], where the axial vibrations to the first order have been considered, too, and the coupling between axial and transversal oscillations is deeply investigated. Interesting, and complex, behaviour charts have been reported, showing that coupled and uncoupled solutions can coexist or not, both being hardening or softening.
The difference obtained by considering alternative definitions of the curvature (geometric, i.e. dθ dS , vs mechanical, i.e. dθ dZ , dS and dZ being the infinitesimal elements of the deformed and undeformed configurations, respectively) in the constitutive behaviour has been the subject of [9] and [10]. It has been shown that the difference is of the order of some percents for thick beams and is totally negligible for slender beams. While in the previous authors' work free vibrations have been considered, in [11] and [12] the forced vibrations have been investigated. In [11], the multiple time scales method has been used, and again the attention was focused on the dependence of the hardening/softening response on the spring stiffness κ. Analytical results have been compared with numerical (FEM) results, showing a good agreement up to moderately large displacements, according to the fact that analytical solution is valid only up to the third order. The appearance of some superharmonic and internal resonances has been observed in the numerical simulations. Still in the forced regime, the comparison between numerical and experimental results was instead the goal of [12], where a very good matching has been observed.
The frequency response curves of higher order resonances have been discussed in [13], where, among other, it has been stressed that the nonlinear behaviour may strongly depend on the mode order, for example the first mode can be hardening and the second softening.
Within this wide body of research, one aspect has not yet been investigated, namely the internal resonance between axial and transversal modes. Indeed, in [8] it was assumed that axial and transversal modes were far away from internal resonance, so that the axial-transversal coupling was only due to the nonlinear effects occurring in kinematically non-condensed structural models (e.g. [14]). Here, on the contrary, this hypothesis is relaxed, and the internal resonance is fully investigated, both in the linear and nonlinear regimes. In addition to its theoretical and practical effects, this also permits to explain a singularity of the nonlinear correction coefficient observed in [5] and in [11]. As discussed in more detail in Sect. 2.1 forward, a similar singular behaviour of the nonlinear correction coefficient, associated with the transition from hardening to softening and due to an internal resonance, has been earlier highlighted in shallow cables [15,16] and in shells [17,18], where only transversal dynamics have been studied.
A transition from hardening to softening, due to 1 : 2 internal resonance and not necessarily related to a singular behaviour, has been observed also in a pipe con-veying fluid [19], where the varying parameter was the fluid velocity, and in a double beam with a tip coupling spring and magnetic interaction [20], where the driving parameters were the amplitude of the excitation and the damping.
Internal resonance, or modal interactions [21], have been deeply investigated in the past [22][23][24], since this event can be dangerous (if not properly detected [25]) or useful (if properly exploited, for example in the field of energy harvesting [26]). However, it seems that few attention has been paid to the internal resonance between longitudinal and flexural modes. Axialtransversal internal resonances of cables have been investigated, for example, in [27,28], of moving belts in [29], while for beams this phenomenon has been investigated in [30], again in the field of axially moving systems. Very few investigations of internal longitudinaltransversal resonance seem to be available for nonmoving beams [31].
The paper is organized as follows. The mechanical model is illustrated in Sect. 2, including a summary of previous results (Sect. 2.1) where the singularity of the nonlinear frequency ω 2 is highlighted. The occurrence of various axial-transversal internal resonances in the slenderness-spring stiffness parameters space is investigated in Sect. 3 in the linear regime. Then (Sect. 4), attention is focused on 1:2 internal resonance, but now a nonlinear dynamic analysis is developed with the multiple time scales method. The main results are reported in Sect. 5, while a comparison of analytical and numerical outcomes is presented in Sect. 6. The paper ends with some conclusions and suggestions for further developments (Sect. 7).

The mathematical model
We consider the Timoshenko beam illustrated in Fig. 1. Its rest configuration is rectilinear, along the Z direction, it is made of a linearly elastic and homogeneous material, and its cross section is constant. W (Z , T ), U (Z , T ) and θ(Z , T ) are the axial and transversal displacements of the beam axis, and the rotation of the cross section. Z is the spatial coordinate in the reference (straight) configuration, measuring the physical distance from the left boundary, and T is the time.
The axial, shear and bending stiffnesses are E A, G A and E J , respectively. The axial and transversal mass per unit length is ρ A, while the rotational inertia is ρ J . κ is the stiffness of the linear spring at the right-end of the beam (Fig. 1), which acts in the Z (axial) direction. C W , C U and C θ are the (linear) damping coefficients along the respective directions. The beam is excited by a dead load P U (Z , T ) = F(Z , T ) acting in the X direction, which is perpendicular to Z . Introducing axial, P W (Z , T ), and rotational, P θ (Z , T ), loads is conceptually easy but will not be pursued to limit the computational efforts.
The kinematically exact equations of motion have been derived in [6] (see also [5,8,11]) and are given by: where the prime denotes derivative with respect to Z and dot derivative with respect to T . The boundary conditions in the transversal direction are: where M is the bending moment, i.e. the geometrical curvature is considered. A possible alternative is to use the mechanical curvature, assuming M = E J dθ dZ [9,10]. The boundary conditions in the axial direction are: where H o (Z , T ) is the internal horizontal force in the Z direction and is given by

Previous results
By applying the Poincaré-Lindstedt method (in [5,6]) and the multiple time scales method (in [11]), and by considering only the transversal displacements up to the first order, the following approximate solution was previously obtained for the free oscillations: where U a is the amplitude of the transversal motion, n ∈ N the mode number, and where, most importantly, is the nonlinear frequency of the free motion. In (7), ω 0 is the dimensionless natural (linear) frequency of the transversal motion, while ω 2 is the dimensionless nonlinear correction coefficient, also known as "nonlinear frequency correction", which measures how the nonlinear frequency is affected by the amplitude of the motion. It was the most important result, since it summarizes the nonlinear behaviour of the system: hardening for ω 2 > 0, softening for ω 2 < 0.
For the boundary conditions (2), it is possible to compute where n ∈ N is the order of the transversal (flexural) natural frequency and the following three dimensionless parameters are introduced: (l is the slenderness of the beam, ν the Poisson coefficient and χ the shear correction factor) so that κ h is the dimensionless stiffness of the axial spring at Z = L, that will be used in the rest of the paper. We consider ν = 0.3 and χ = 1.17, namely z = 0.3287, so that the unique parameters are l and κ h , together with the order n of the natural frequency. For the same boundary conditions, the nonlinear correction coefficient is where the parameters appearing in (11) are reported in "Appendix".
For κ h = 50 and n = 1, the function ω 2 (l) is plotted in Fig. 2, which corresponds to Fig. 10b of [5]. For the scope of the present work, the most relevant property of Fig. 2 is that ω 2 (l) has a singular point for l sing = 8.149112. It has also a zero point at l zero = 11.379729. When crossing these points, the sign of ω 2 changes, and, for increasing l, we pass from hardening to softening and again to hardening.
The singular and zero points are not specific of the considered value of κ h , but they persist, as shown in Fig. 3, corresponding to Fig. 12c of [5].
In [5], eq. (27), looking at the expressions of ω 2 , it was noted that the singular points occur when namely when the denominator of c 2 is zero. This phenomenon was also observed in [11], and in [32] it is mentioned that it is due to an internal resonance between axial and transversal modes. Confirming this property in the considered mechanical system, and studying in depth the internal resonances, is the main goal of this paper. Before to proceed, we note that the singular points in Fig. 3 occur for quite low values of the slenderness (thick beams, but still compatible with the Timoshenko beam theory). Actually, for n = 2 the singular points happen for larger values of l (in the realm of slender beams), as shown by the example of Fig. 4, which corresponds to Fig. 3 of [11] (but note that a different For n = 2, the zero and singular points are reported in Fig. 5. Note that in this parameters window there are two singularities and zeros. The lowest zero and sin- Although the motivation of the present paper comes mainly from the authors' previous work summarized above, it is worth mentioning that the singular behaviour of the nonlinear correction coefficient in correspondence of internal resonances, and in particular 1 : 2 internal resonance, has been previously observed for different mechanical systems in the framework of the reduction methods evaluation. In [15], the singularity has been reported for the first in-plane (vertical) mode of a suspended cable. Here, the varying parameter is an elasto-geometric parameter taking into account stiffness and sag-to-span ratio, and the singularity is due to the breakdown of modal expansions not explicitly considering the internal resonance. For the same structure, a similar, but deeper, analysis has been reported in [16], where also the out-of-plane behaviour has been considered. Here, it is shown that discretized approaches, where the Galerkin reduction method is applied before the MTS, fail to detect the singular behaviour if a sufficiently large number of modes is not considered.
Using nonlinear normal modes to highlight hardening/softening behaviour of nonlinear systems [33], the singular behaviour of the nonlinear correction coefficient for a free-edge shallow spherical shell has been shown in [17], where the driving parameter is the aspect ratio, which is governed by the geometrical properties of the system. Within a study of the effect of the number of retained modes in a Galerkin reduction, many curves similar to Figs. 2 and 4 have been reported. It is remarked that only 1 : 2 internal resonance of certain modes are able to generate the singular behaviour, while other modes and other resonances (e.g. 1 : 3) do not have this effect. In [18], it has been shown that the addition of damping will smooth the singular behaviour, while keeping the large (but finite) values of ω 2 for small values of damping.
In the previous papers, the singular behaviour, and the 1 : 2 internal resonance lurking in the background, have been reported. However, a detailed asymptotic expansion considering the internal resonance is not carried out. This analysis has been done in [19,20], where, however, no singular behaviour is observed in the background. In the present work, we complete, and somehow connect, the two groups of works, while referring to the longitudinal-transverse coupling in the mechanical system of our interest.

Longitudinal-transversal internal resonance: linear analysis
The occurrence of the longitudinal-transversal internal resonance in the parameters space (κ h,l ) can be detected in the linear regime, where the two modal problems are independent from each other (since the beam is rectilinear). This constitutes the aim of this section. The nonlinear coupling in the neighbourhood of internal resonances will be investigated in the next section. The transversal natural frequency ω f 0 is given (in dimensionless form) by (8). The axial natural frequency ω a0 , the solution of can be easily computed by standard methods (see also Sect. 4.1) and is given by x being the solution of The associated mode shape is w(Z ) = sin(x Z/L). Note that for very small values of κ h /l 2 the solutions of (15) are x m ∼ = π/2+(m−1)π , m ∈ N being the order of the axial natural frequency, while for very large values of κ h /l 2 the solutions of (15) are x m ∼ = mπ . Thus, for varying κ h /l 2 ∈ [0, ∞[ the following bounds hold: The expression (15) can be rewritten as (16) is identical to (12), and this clearly shows that the singular points of ω 2 are due to the 1 : 2 internal resonance between transversal and axial modes. Reasoning in a similar fashion, it is not difficult to foresee that singular points of ω 4 (the next term in the development (7)) will be due to a 1:4 internal resonance between transversal and axial modes and so on. This aspect will not be pursued here and is left for future works.
To generalize the previous considerations, we investigate the 1 : ir internal resonances, those such that The study of other internal resonances, for example 2 ω a0 = ω f 0 , is left for future works. Combining (7), (8), (14), (16) and (17), it is easy to obtain For ir = 1; 2; 3, the loci of internal resonance in the parameter space (κ h , l) are reported in Fig. 6, where n is the order of the transversal natural frequency, see (8). To understand the order m of the axial frequency, for each considered value of l, the lowest value of κ h is associated with the first axial mode, the next to the second and so on. For example, for ir = 1 (Fig. 6a) and for l = 20 we have that for κ h = 75.0464380 (point A) there is a 1:1 internal resonance between the second (because n = 2) transversal mode and the first axial mode (so that m = 1); for κ h = 829.704872 (point B) there is 1:1 internal resonance between the fourth (because n = 4) transversal mode and the second axial mode (so that m = 2).
In Fig. 6b, the point C corresponds to the case of Fig.  2: here, we have a 1 : 2 internal resonance between the first (n = 1) transversal mode and the first axial mode (m = 1). The point D corresponds to the case of Fig.  4: here, we have a 1 : 2 internal resonance between the second (n = 2) transversal mode and the first axial mode (m = 1). Note that, for the same value of κ h , there is another 1 : 2 internal resonance between the second transversal mode and the first axial mode on the lower n = 2 curve, for l = 5.7977.
We notice that internal resonances are more sensitive to the slenderness l than to the spring stiffness κ h . In fact, at least in the considered range of parameters and irrespective of the value of ir, internal resonances occur only for "fixed" values of l (or in narrow intervals of l), and there are large intervals of l (sort of "band gaps") without internal resonances. On the contrary, they occur for every value of κ h .
To end this section, we note that by varying the transversal boundary conditions (for example assuming fixed instead of hinged left boundary condition, (2)), the axial frequencies do not change, while the transversal ones are shifted. This permits to conclude that the same internal resonances observed in this section are expected to exist also for other boundary conditions, and they are only translated in the parameters space. The check of this guess is left for future work.

1 : 2 internal resonance: nonlinear analysis
In this section, we study, in the nonlinear regime, the 1 : 2 internal resonance illustrated in Fig. 6b, via an asymptotic analysis up to second order. In addition to its interest in explaining the singularities of the nonlinear correction coefficient ω 2 , it has an interest "per se", since this resonance has interesting dynamical behaviours, as we will see. The detailed nonlinear investigations of other internal resonances can be done analogously.
The multiple time scales (MTS) method [34] is used, by considering to the first order both axial and transversal modes. Also, in [11] the MTS was used, but there only the transversal mode is considered to the first order, so precluding the possibility to detect internal resonance (that was not the goal of that paper, indeed). Axial and transversal modes to the first order are also considered in [8], but there it is implicity assumed that they are not in internal resonance, since the aim was to detect the axial-transversal nonlinear coupling in a general fashion, and not only around specific values of the parameters where internal resonance occurs.
With the MTS method, the solution is sought after in the asymptotic form where ε is a bookeeping small parameter introduced to stress that we are studying moderate displacements and rotation around the rest position. T 0 = T is the physical (fast) time, and T i = ε i T , i ≥ 1 are the slow times.
It is further assumed that the excitation is harmonic in time (with frequency ), and that damping and load scale according to Inserting the expressions (19)- (20) in the asymptotic expansion of the governing equations (1), and equating to zero the coefficients of ε i , a sequence of linear problems is derived. They are investigated in the next subsections.

First-order problem
The first-order equations read and the related boundary conditions are (the dependence on time is omitted for simplicity) The solution of (21), (22) is given by (the over bar stands for complex conjugate and I is the imaginary unit) where the functions W 1a (Z ), U 1a (Z ) and θ 1a (Z ) are reported in "Appendix", and A(T 1 ) and B(T 1 ) are complex amplitudes of the transversal and axial oscillations, respectively. ω a0 is given by (14) (with x given by (15)), while ω f 0 is given by (7) with ω 0 given by (8). Note that in [5][6][7]9,11], but not in [8], it was assumed W 1a (Z ) = 0.

Second-order problem
The second-order equations read (24) and the related boundary conditions are (the dependence on time is omitted for simplicity) It is now necessary to introduce the 1 : 2 internal resonance condition, together with the external resonance, where σ i and σ e are the detuning parameters measuring the distance (in the frequency) from the perfect internal and external resonances, respectively. Note that, combining (26) and (27), we have so that the detuning between 2 and ω a0 is given by The particular solution of (24) is given by It is worth to note that the homogenous solution of (24), that adds to (29), is not needed in this work. Inserting (29) in (24), and in the boundary conditions, we obtain the following problems for W 2a (Z ), U 2a (Z ) and θ 2a (Z ): where the functions f W , f U and f θ are reported in "Appendix".

Frequency response curves
The solvability conditions for the second-order problems (30) are where the parameters r 1 , . . . , r 6 are reported in "Appendix". Note that the load is accounted for in r 6 . Remembering that x = x(κ h , l, m), we observe that the other coefficients depend on l, κ h , z (that is kept fixed in this work) and the order of transversal (n) and axial (m) natural frequencies.
As customary [34], the solution of (31) is sought after in the polar form where a(T 1 ) and b(T 1 ) are real value amplitudes of the transversal and axial (first order) motions, respectively, and β a (T 1 ) and β b (T 1 ) are real value phase differences. In fact, rearranging (23) we obtain Inserting (32) in (31) and separating real from imaginary parts, we obtain, after some algebra, that are commonly known as modulation equations.
Here, the dependence on T 1 is omitted for brevity. Because of (33), steady oscillations of the beam correspond to equilibrium points of (34). Setting equal to zero the derivatives in (34), we obtain an algebraic system in the four unknowns a, b, β a and β b , that are now constant.
From the third and fourth equations, we obtain b = r 2 2ω a0 r 1 Using (35) in the first and second equations of (34), we obtain tan(β a ) =

4(c
and the equation which is third order in the unknown y = a 2 , so that its closed form solution is known. Once a(l, κ h , n, m, c W , c U , c θ , r 6 , ) has been determined, b, β a and β b can be computed from (35) and (36). When all parameters are fixed, and only the external excitation is varied, the searched frequency response curves (FRCs) are obtained.
Equation (37) is cubic in a 2 . Thus, there always exists a real solution. Furthermore, it may happen that three real solutions exist. This occurs in particular in the neighbourhood of resonance, as we will see in the following.
In the very special case of perfect internal (σ i = 0) and external (σ e = 0) resonances, i.e. ω a0 = 2ω f 0 = 2 , Eq. (37) reduces to a 1 16 In this case, the asymptotic development of the solution with respect to c W is given by which shows that for c W → 0 we have b → ∞. This highlights the role played by c W in the resonance. The stability of the solutions previously obtained can be determined by computing the eigenvalues of the Jacobian of the right-hand side of (34) at each equilibrium point. If all the four eigenvalues have negative real part, the solution is stable; otherwise, it is unstable. In particular, if there is a real positive eigenvalue, the solution is a saddle, while if there is a complex eigenvalue with positive real part, the solution is a source.
By varying the parameters, if a real eigenvalue passes from negative to positive, we have a saddlenode bifurcation, while if a complex eigenvalue passes from negative real part to positive, we have a Neimark-Sacker, or secondary Hopf, bifurcation, and a quasiperiodic solution is born in the mechanical system. Because of their geometrical properties, saddle-node bifurcations can be detected by looking for the zeros of the discriminant of Eq. (37), where we pass from 1 to 3 real solutions.

Backbone curves
In the absence of excitation, r 6 = 0, and damping, c W = c U = c θ = 0, the beam undergoes a free motion, of frequency (see 33). In this case, (37) simplifies to which can be rewritten as Thus, excluding the trivial solution a = 0, we have and then (the sign function is sign They represent the analytical expressions for the socalled BackBone Curves (BBCs) of coupled oscillators, which give the amplitudes of the oscillation as a function of the vibration frequency . As it happens for the BBCs of uncoupled oscillators [34], they are also the loci of maximum points of frequency response curves, as we will see in the next section.
When all parameters but are fixed, only σ e is varying in (42) and (43), so that a(σ e ) ∼ √ σ e (2 σ e − σ i ) and b(σ e ) ∼ σ e sign(2 σ e − σ i ), this being a property that helps to understand the behaviour of the BBCs. In particular, we note that a(σ e ) exists only when σ e / ∈ [min{0, σ i /2}; max{0, σ i /2}], and locally it behaves like a square root. b(σ e ) exists in the same range, where it is (piecewise) linear.
In the very special case of perfect resonance, σ i = 0 and ω a0 = 2ω f 0 , the previous expressions simplify to and both curves become proportional to |σ e | (see forthcoming Fig. 9d).

Results
Although the analysis of the solution can be done in dimensionless form, to refer to a real case and to compare the results with those obtained by numerical simulations, we prefer to deal with a dimensional case. We choose a beam of 0.05 × 0.05 m 2 square cross section, made of steel (E = 2.1×10 11 N/m 2 , ρ = 7850 kg/m 3 , ν = 0.3), of length L = 0.5 m. We then have l = 10 √ 2 = 34.6410 and z = 0.3287, so that we are exactly in the case of Fig. 4. We also have ω f 0 = 11089.77 rad/s (the period is 5.66 × 10 −4 sec) and ω a0 = 10344.39 x rad/s. For this case, the perfect 1 : 2 internal resonance between the second transversal mode (n = 2) and the first axial mode (m = 1) is obtained at κ h,ir = 1661.24, see point "D" in Fig. 6b. We will vary κ h around this value.
We further assume C W = C U = 50 N sec/m 2 and C θ = 1.5 N sec. The load is a concentrated force Q applied at L/4 (in order to excite the second transversal mode), so that r 6 = Q sin(n π/4)/2 = Q/2 N (since n = 2). Any other distribution of load providing the same r 6 is equivalent to the present one. For comparison with the forthcoming dynamical case, we note that the transversal static displacement in the point where the force is applied is, within the small displacements linear theory, U (L/4) = (3/256)(Q L 3 /E J ) = 0.0134 mm for Q = 1000 N.
Note that, according to (14) and (15), ω a0 depends on κ h , and thus, it slightly varies when κ h varies from the value 1620. The FRCs and the BBCs are reported in Fig.  7, which shows how the excitation amplitude increases the amplitude of the response, as expected, but not in a linear way. The dynamic transversal displacement in correspondence of the peak is about 100 times larger than the static one. The stability of each branch is determined as illustrated at the end of Sect. 4.3. For example, for Q = 1000 N and = 11120 rad/s we have the following equilibrium points and associated eigevalues: Overall there are two unstable solution branches of saddle-type (reported in red), that are bounded by two saddle-node (SN) bifurcations each, and an unstable solution branch, of source-type (reported in blue), that is bounded by two Neimark-Sacker (NS) bifurcations, that are seen in Fig. 7. For Q = 1000 N, the maximum real part of the eigenvalues is reported in Fig. 8 for each branch, from which we clearly see the stable (max Re(λ) < 0) and the unstable (max Re(λ) > 0) solution branches. The bifurcation points are where these curves cross the abscissa axes ((max Re(λ) = 0)), that happens for The source-type unstable solution branch is quite narrow. For low values of the load (Q = 300 N), there is a unique stable solution branch, having peaks in correspondence of the resonances.

The transition from hardening to softening
The transition from hardening to softening is illustrated in Fig. 9 to show what really happens through the (apparent) singularity of Fig. 4. For κ h quite below κ h,ir (Fig. 9a), we are far from the internal resonance. In correspondence of ω a0 /2, there is no amplification on U , because the excitation is not in primary resonance with the transversal mode, while we have a classical resonance curve in W , which is practically vertical.
Here, the resonance mechanism is the following: (1) the external excitation induces small vibration on the transversal displacement; (2) the weak nonlinear coupling from U to W , illustrated in [8], induces small oscillations on the axial displacement; (3) since the external excitation is just in order-2 superharmonic resonance with the axial mode, there is a resonance amplification, linear because the transversal mode transfers only a small amount of energy to the axial one. In turn, around ω f 0 , the FRC of U is hardening (in agreement with ω 2 > 0 in Fig. 4), and that of W behaves similarly, due again to the mere nonlinear coupling away from internal resonance [8]. Note that the transversal displacement (directly excited) is about 15 times larger than the axial one. Increasing κ h (Fig. 9b), the major change is that a softening resonance curve is observed in both W and U around ω a0 /2. The meaningful nonlinear effect of the internal resonance is apparent in the left peaks, the W one being much more important than the corresponding hardening peak around ω f 0 which represents the nearly unchanged effect from U to W due to the nonlinear coupling. Note also the clearly visible piecewise linear behaviour of the BBCs of b in Fig. 9b2.
By further increasing κ h (Fig. 9c), thus getting closer to the exact internal resonance, the softening FRC around ω a0 /2 becomes more important, both in terms of U and W . In W , a third peak appears exactly at ω a0 /2, suggesting a further second-order couplinginduced effect from U to W ; however, this peak belongs to the unstable solution branch. The hardening branch of the FRC around ω f 0 is still practically unchanged.
For κ h = 1660 (Fig. 9d), i.e. practically in the perfect resonance, the left softening and the right hardening cusps become symmetrical with respect to the line = ω f 0 = ω a0 /2, and thus, they reach the same importance. The third, unstable, W -peak is clearly evident in Fig. 9d2. Here, the coupling due to the internal resonance is maximum, and the two BBCs, one hardening and one softening, ensue from ω f 0 , and are piecewise linear, as found in (44).
The increment of κ h (Fig. 9e) implies that the FRC branch associated with ω f 0 becomes the left one, which is softening, while its hardening branch (associated with ω a0 /2) becomes slightly less important in terms of U , but more important in terms of W . Indeed, note that Fig. 9e is "specular" to Fig. 9c. This trend is confirmed by Fig. 9f, which is "specular" to Fig. 9d, and by Fig. 9g, which is "specular" to Fig. 9a.
The conclusion is that, taking properly into account the internal resonance, there is no singularity in the hardening to softening transition by varying κ h . There is, instead, a sudden (discontinuous) jump, of finite magnitude, from hardening to softening, and in the discontinuity point κ h,ir , where exact internal resonance occurs, the BBCs are linear instead of proportional to a square root.
Finally, we report in Fig. 9h  show how the softening behaviour decreases (the U curve becomes almost vertical, and thus, the nonlinear correction coefficient almost zero). This is in agreement with the behaviour illustrated in Fig. 4, and it is somehow surprising, since the present analysis is expected to be reliable only in a rather strict neighbourhood of κ h,ir .

An isolated branch of solution
The born of the left resonance curve (see Fig. 9b) is actually more involved than expected, and interesting. An isolated branch of solution, also termed as "isola", is born at about κ h = 1580 (Fig. 10a, b). By increasing κ h , this closed branch enlarges (Fig. 10c, d) and finally touches, at about κ h = 1591.5 (Fig. 10d, e), the nonresonant, small amplitude, branch. Then, a classical resonant (softening) FRC is observed (Fig. 10f). The same phenomenon has been observed, for example, in [35] (see their Fig. 6) and in [36] (see their Fig. 2). It is just for isolated branches that analytical methods, even approximate, are necessary: here multiple time scales are used, but in other papers (e.g., [37]) isolas have been detected by the harmonic balance method. In fact, they cannot be obtained by a path following algorithm starting from the main solution branch unless varying another parameter (for example, the stiffness of the spring), besides the excitation frequency, which may allow to reach isolas in an extended parameter space. If varying only the excitation frequency, they can be found numerically only via a suitable and lucky choice of initial values in either the original governing equations or the reduced (e.g., modulated) ones, and then possibly characterized systematically by building basins of attraction, which is however difficult and demanding, especially for high dimensional systems. Detecting stable isolated solutions is not only important from a theoretical point of view. In fact, even if these attractors may have a small basin of attraction, and thus be considered "minor" or "rare", it may happen that their basins strongly erode, by fractality, the basin of the main attractor (the one used in applications), that thus loses robustness and reduces its safety, leading to instability for moderately large perturbations, even if the attractor is stable in the Lyapunov sense [38].

Numerical simulations
To confirm the previous analytical findings, finite elements method (FEM) simulations are performed. The same beam illustrated in Sect. 5 is considered.
The commercial software Abaqus_CAE © is used. The beam is discretized using 100 B21-type elements. An explicit direct time integration is used, with time step T = (2π)/( × 80) sec and, to give up the transient vibration, each excitation lasts 2500 periods. Furthermore, in the solver a double precision is used together with full nodal output precision.
Relatively small changes in (1 rad/s) are used in the frequency sweeping (both forward and backward), and the solution (shape deformation and rotation angle) obtained in the previous step is used as initial condition, as done in [11].
Axial and transversal dampings proportional to the mass matrices are considered, with the same values of the damping coefficients used in Sect. 5.
To further confirm our results, we have done computations also with another commercial software, Midas_GEN ©, but since the outcomes are comparable we do not report these simulations.
First, the axial and transversal natural frequencies have been computed in the linear regime, verifying that they perfectly match those computed analytically. Then, the FRCs have been obtained in the geometri-cally nonlinear regime, i.e. with large displacements. Since a brute force following algorithm has been used to determine the numerical FRC, for increasing and decreasing frequencies, only the stable branches have been detected. To draw also the unstable branches, a continuation algorithm should have been used, which is, however, out of the scope of the present work.
The comparison between analytical and numerical results is reported in Fig. 11, for a fixed value of κ and for increasing values of Q. We conclude that an excellent agreement is obtained on the main (i.e. primary resonance) stable solution branches. Surprisingly, the peaks of the numerical curves are higher than those of the analytical ones, this being likely due to a different fine tuning of the damping, that has large effects on the peak value. Note that no tuning of damping is done, and the nominal values are considered in numerical and analytical simulations.
Around ω a0 /2, we observe that the numerical solution is not able to catch the upper, resonant, stable branch. This is due to the source-type unstable solution branch that breaks the main resonant solution branch and leads to a jump of the theoretical solution. The numerical jump is not observed, likely because the solutions have very small basins of attraction.
Keeping Q = 1000 N fixed, and varying the spring stiffness, the comparison between analytical and numerical results is reported in Fig. 12. In addition to the overall excellent agreement, with some minor quantitative differences on W (which can be reduced by tuning the damping), the most important result is that through the Neimark-Sacker bifurcation a quasiperiodic attractor is born, according to what the theory predicts for this event [39]. It is highlighted by the cloud of points obtained for each frequency by sampling the maximum amplitudes of the last 500 periods of the excitation, to give rid of the transient behaviour.
This is a further confirmation of the reliability of the proposed analysis.

Conclusions and further developments
The internal resonance between axial and bending modes in a Timoshenko beam with a boundary axial spring has been investigated, with the aim (1) of studying the ensuing modal interaction in a further mechanical system, (2) of explaining a singular behaviour previously reported in the literature, and (3)   Comparison of numerical (circles) and analytical (continuous lines) FRCs for κ h = 1620 (κ = 1.4175 × 10 9 N/m). In this case, ω a0 /2 = 11041.10 rad/s is slightly lesser than ω f 0 = 11089.77 rad/s. a, b Q = 800 N; c, d Q = 1000 N; e, f Q = 1200 N the transition from softening to hardening through the 1 : 2 internal resonance by varying the axial spring stiffness.
It has been shown that internal resonance, while structurally unstable, is a robust phenomenon, since it occurs for all values of the end stiffness, and for "any" value of the resonance ratio ir. It is only needed to properly select the beam length.
Attention has then been paid to the nonlinear behaviour around the 1 : 2 internal resonance. The fre- 11089.77 rad/s; c, d κ h = 1659.4 (κ = 1.4520 × 10 9 N/m), so that ω a0 /2 = 11087.65 rad/s and ω f 0 = 11089.77 rad/s, namely we are practically in exact internal resonance quency response curves, and the associated backbone curves, have been obtained analytically by the multiple time scales method. It has been shown that, by varying the spring stiffness across the resonance value, the frequency response curves of the internally resonant axial and transversal modes intersect and exchange with each other, leading to an instantaneous, but not singular, jump from hardening to softening of the flexural mode, thus explaining an apparent singularity previously reported when internal resonance is not taken into account.
A detaching of an isolated solution branch from the frequency response curves, similar to others previously reported in the literature, has been observed, too.
Finally, the analytical frequency response curves have been compared with their counterparts obtained by finite element numerical simulations. Overall, a very good agreement has been observed for the considered values of the parameters.
As far as further developments are concerned, at least the following ones can be mentioned: • To explain possible singular behaviour of higher order nonlinear correction coefficient ω 4 in terms of 1:4 internal resonance; • To study other internal resonances, including a combination one, between axial and transversal modes, both in linear and nonlinear regimes; • To consider different boundary conditions; • To exploit the investigated internal resonance for practical purposes (e.g. energy harvesting) and engineering design; • To perform experimental investigations tailored around internal resonance; • To demonstrate that the coupled backbones tend to their uncoupled counterparts when getting far from the resonance.