Nonlinear dynamics of a tubular beam considering distortion of the cross sections and internal resonances

The nonlinear dynamics of a thin tube under the action of a harmonic external load is addressed in the paper. Use is made of a beam-like model which extends the Timoshenko beam model with further kinematic descriptors, related to the change in shape of the cross section. The external load is applied in half cap of the pipe, directly triggering both the bending of the axis-line and the flattening of the cross sections. The equations of motion are projected on a reduced basis constituted by the first three linear modes, and then the solutions are sought via the multiple scale method, for two different external resonance conditions. Internal resonances among the modes are considered as well. The outcomes, compared with pure numerical solutions, highlight the possible energy exchange between local modes, i.e., those describing flattening and warping of the cross sections, and global modes, i.e., those related to bending of the axis-line and rotation of the cross section of the pipe.


Introduction
The use of tubes as structural members is widespread in lots of engineering applications, ranging from civil to industrial, aerospace and many other contexts. Hence, the evaluation of their carrying capacity is a compelling step in the design development, as the typical thin-walled nature of pipes makes this aspect a key point. Indeed, demanding cares are requested to consistently deal with the possible lack of validity of the Saint Venant principle and the significant contribution of the distortion of the cross sections. In this framework, the possibility to model tubes as beams or beamlike structures would represent an undoubted asset, as compared to more demanding bi-dimensional or threedimensional continuum theories. However, classical beam theories like Euler-Bernoulli and Timoshenko [1] require to be enriched, in order to overtake the hypothesis of rigid cross section.
For instance, the Vlasov theory [2], introduced to address non-uniform torsion of tubular beams, serves as main instrument in including the effect of warping of cross sections in shaft models, as well as in providing advanced and reliable contributions in the mechanics of pipes. As a further example, the nonlinear interaction between bending of pipes and flattening of their cross sections is the main focus of the Brazier theory [3,4], which provides physical explanations and functional tools to engineers to address the softening behavior of bent tubular beams. In some cases, soft elastic cores are included to prevent flattening [5].
Many other efforts have been made in the last decades by scholars in developing high-level beam theories [6][7][8][9]. The generalized beam theory (GBT) lies in this line of research [10][11][12][13][14], introducing linear combinations of assumed shape functions to describe bending, torsion and cross-section distortion of thin-walled beams. Recently, in [15], the Euler-Bernoulli beam model was endowed with descriptors for the distortion of the cross section, to deal with multi-layered pipes under flexural static actions. Lately, in [16,17], the same idea, which originally comes from GBT, was broadened to the Timoshenko beam model, including effects of flattening and warping of the cross sections as assumed shapes, amplified by unknown variables. The same pipe model was then used to analyze the nonlinear dynamic response in [18], after consistent evaluation of the inertial contributions. In that specific case, the nonlinear coupling came both from stiffness and inertial terms, and triggered internal resonances between modes, which were related both to global (bending) and local (cross-section distortion) behavior. External resonance due to a load covering the whole cross section was considered as well, and the effects of the softening contribution provided by the cross-section change in shape were analyzed.
In this paper, starting from the pipe model proposed in [18], a different external load is considered: Here, it is assumed to be distributed only on half cap of the cross sections. This specific aspect induces a direct loading toward the flattening modes of the cross sections of the tube. Hence, 1:1 external resonance conditions with one of the local modes, combined with internal resonances between local and global modes, may potentially cause energy transfers from cross-section distortion to bending, occurrence which is worthy of investigation. The analysis is addressed by the multiple scale method (MSM), applied to the equations of motion which are made discrete by a Galerkin projection. Two different implementations of the MSM are carried out, depending on the local mode involved in the considered external resonance condition. Numerical integration of the equations of motion are used to compare and validate the asymptotic solutions.
The paper is organized as follows. In Sect. 2, the beam-like model is briefly described, in Sect. 3, the discretized nonlinear equations of motion are obtained via a Galerkin projection, in Sect. 4, the two different implementations of the MSM are described, and in Sect. 5, the numerical results are presented and discussed. Finally, the conclusions are drawn in Sect. 6.

Model description
The formulation of the beam-like model used here to address the nonlinear dynamics of a pipe with thin annular cross section is extensively described in [18]. For the sake of completeness, here its main features are only briefly recalled, leaving the details in [18], but highlighting the differences related to the load and resonance conditions. An in-plane Timoshenko beam model is introduced (Fig. 1a), as constituted by a straight axis-line spanned by the abscissa s ∈ [0, l] in directionā 1 , where l is the initial length, and by infinite initially transverse cross sections, parallel to the plane spanned byā 2 ,ā 3 . The unitary vectorsā 1 ,ā 2 ,ā 3 are mutually orthogonal. The beam is clamped at s = 0 (cross-section A) and free at s = l (cross-section B). The kinematic variables are u(s), v(s), ϑ(s), where the first two variables representā 1 -andā 2 -components of the displacement u of the axis-line points, respectively, whereas the third one describes the cross-section rotation aboutā 3 . Moreover, as an extension of the Timoshenko beam model, further kinematic variables are introduced, referred to as a p (s), a w (s), which describe in-plane and out-of-plane change in shape of the cross section, respectively. The physical meaning of a p (s), a w (s) comes from an identification procedure of the beam-like model through a three-dimensional realization of the pipe, seen as an assembly of infinite longitudinal fibers and transversal ribs (Fig. 1b), and having length l, mid radius R, and thickness h R. In particular, a p (s), a w (s) turn out to be the amplitudes of assumed flattening (Fig. 2a) and warping (Fig. 2b) shapes of the annular cross section of the pipe at the generic abscissa, respectively.
The strain measures of the Timoshenko beam are consistently introduced: the longitudinal strain ε 0 (s), the transversal strain γ 0 (s), the bending curvature κ 0 (s), as well as strain components relevant to the cross-section change in shape, namely α p (s), β p (s), α w (s), β w (s), referred to as local components. Hence, the nonlinear strain-displacement relationship, seriesexpanded up to the third order, is:  Local distortion on the cross section: a assumed shape for the flattening and amplitude a p ; b assumed shape for the warping and amplitude a w where prime stands for s-derivative. Boundary conditions for clamp at cross-section A read: The virtual work theorem allows one to determine the weak form of the dynamic equilibrium equations.
More specifically, the internal virtual work for the beam-like structure reads: where δ is the variational operator, T , M are the shear force and bending moment of the planar Timoshenko beam, D j , B j are distortion and bi-distortion force components, dual to the local strain components, for j = p, w, and λ is a Lagrange multiplier, introduced in order to nullify the longitudinal strain ε 0 , as it is usual in case of cantilevers [19][20][21][22][23][24]. Moreover, the external virtual work for the beam-like structure reads: where f u , f v , c, g j represent the external distributed forces and couples, work-conjugate of the generalized displacements, andf u ,f v ,c,g j are the inertial counterparts. Substitution of Eq. (1) in Eq. (3), imposition of the virtual work equation δW int = δW ext , for all kinematically consistent δu, δv, δϑ, δa p , δa w , δλ, provides the weak form of the dynamical equilibrium equations.
The constitutive law in case of linear elastic material of Young modulus E and transversal elastic modulus G is obtained after the application of the identification procedure from the three-dimensional model and assumes the following expression: where the elastic coefficients are: Here, the Poisson ratio is assumed as ν = 0, in order to highlight the pure effect of the coupling between bending and flattening. It is worth mentioning that the three-dimensional model used for the determination of the constitutive law (5), sketched in Fig. 1b, is assumed to allow extension and shear deformation of the longitudinal fibers, as well as bending of the transversal annular ribs. More details on this aspect are given in [16,18].
Consistently, the expressions for the inertial forces and couples are identified as well: with the coefficients: π hρ R, m 5 = 5 8 π hρ R, m 6 = π hρ R, where the dot stands for differentiation with respect to time, indicated as t.
In Eq. (7), the condensation of the variable u is applied, as a consequence of the condition ε 0 = 0 which provides, by Eq. (1-1): Correspondingly, the expression of the Lagrangian multiplier is obtained as well: The identification procedure also provides the expression for the external forces. Here, a time-dependent load per unit volume b v = b 0 cos(Ωt), in the direction a 2 , is uniformly applied on the upper cap of the cross sections, as shown in Fig. 3. Therefore, the load condition is different than the one applied in [18], where the load was applied in the whole cross section. Hence, in the analyzed case and with reference to Eq. (4), the load provides both f v and g p components in the beam-like model, of expression: f v = f 0 cos(Ωt), g p = − 4 3π f 0 cos(Ωt), with f 0 = π h Rb 0 whereas f u = 0, c = 0, g w = 0. In other words, the load produces nonzero work both in the transversal displacement and in the flattening component of displacement.
If free linear vibrations are sought, Eqs. (59)-(61) are written retaining only linear terms and neglecting the external forcing contributions, namely: and at B: It is worth noticing that Eqs. (15)-(23) are uncoupled in the global (v, ϑ) and local (a p , a w ) problems, since coupling only occurs through nonlinear terms. As a consequence, linear modes for the global problem (i.e., those of the Timoshenko beam) are unmodified by the local motion and vice versa. Furthermore, a class of local modes, i.e., involving a p and a w only, is obtained.

Reduced-order model
A Galerkin projection of the nonlinear problem is performed here, using as trial functions the first three modes of the linear problem Eqs. (15)-(23), where one is global (frequency ω 1 ) and two are local (frequencies ω 2 , ω 3 ). Moreover, the frequencies of the higher modes are assumed quite far from the considered ones, so as to neglect their contributions in the response. These assumptions will be lately fulfilled in the numerical example.
The following expressions are hence introduced: ⎛ where q j (t), j = 1, 2, 3 represent the unknown time-dependent amplitudes, and (φ v,1 (s), φ ϑ,1 (s)), (φ p,k (s), φ w,k (s)), k = 2, 3 are the modal components. Substitution of Eq. (24) in the virtual work equation, calculation of the integrals in ds and collection of the terms multiplying δq j , j = 1, 2, 3, produces the reduced ordinary differential equations of motions. In the state space form, they appear as: According to the choice of the trial functions and their normalization, K = diag(ω 2 j ), j = 1, 2, 3 is the (diagonal) stiffness matrix listing on its diagonal the square of the natural frequencies; a linear damping operator C = diag(2ζ j ω j ) is inserted in Eq. (25), being ζ j the damping factors. The load column vector F is defined as: and N is the column vector collecting the quadratic and cubic nonlinear terms: where the single functions are explicitly defined in Appendix B.
variables are expressed as series expansion, after introducing the small scaling parameter 0 < 1: where the different time scales are defined as t 0 = t, t 1 = t, t 2 = 2 t. The linear damping ratios are assumed to be small so that they appear directly at the highest order, i.e., ζ j = 2ζ j (tilde is omitted in the follow).
Internal resonance conditions are considered as well, namely ω 2 2ω 1 , ω 3 3ω 1 in order to possibly address energy exchange between global and local modes. Furthermore, as the external action provides a direct excitation also toward the local modes (Eq. (26)), the following two load cases are considered, inducing different external resonances: To analyze the aforementioned cases, two distinct perturbation schemes are developed to take into account the proper external detunings. It is worth noting that the case Ω ω 1 was addressed in [18].

Case 1: Ω ω 2
For the specified case, the following scaling is adopted for the forcing terms defined in Eq. (26): so that the resonant term will appear at the cubic order, while the non-resonant forcing terms at the linear order. Moreover, external detuning σ is considered as: and the internal detuning parameters ρ 2 , ρ 3 are defined so that: Under these assumptions, the following perturbation equations are obtained: where ∂ j = d/dt j with j = 0, 1, 2. The nonlinear terms are expressed according to the functions defined in Appendix B, while the forcing terms are defined as follows:

Linear-order problem
The solution of the linear-order problem (32), beside the complementary solution, is characterized by the presence of the particular solution related to F 1 . Accordingly, it is defined by the following expression: where (z k , iω k z k ) T is the k-th right eigenvector of the eigenvalue problem given by Eq. (32) made homogeneous. More specifically, it is: Moreover, cc stands for complex conjugate, and the components of the vectors Λ j are: (38)

Quadratic-order problem
After substituting expressions (36) into Eq. (33), the following solvability condition is imposed to eliminate the secular producing terms: with k = 1, 2, 3, being (iω k z k , z k ) T the k-th left eigenvector and R 2 the right-end side of Eq. (33). From Eq. (39), the following amplitude modulation equations are derived: where the coefficients d j ( j = 1, . . . , 6) are defined in Appendix C.
Equation (40), substituted into Eq. (33), allows one to determine the particular solution of the quadratic order problem that can be written in the following form: where NRT represents the non-resonant terms that are not reported here for sake of brevity, while the vectors w j are defined in Appendix D.

Cubic-order problem
After substituting expressions (36) and (41) into Eq. (34), in order to eliminate secular producing terms, the following solvability condition is imposed: with k = 1, 2, 3 and being R 3 the right-end side of Eq. (34). From the latter expression, the following amplitude modulation equations at the third order are derived: The coefficients appearing in Eq. (43) are explicitly defined in Appendix C and their numerical values are shown in Appendix E, with reference to the case study proposed in Sect. 5.1. The reconstructed amplitude modulation equations in the true time t can be written in the form: where the terms ∂ 1 A k and ∂ 2 A k are given in Eqs. (40) and (43), respectively. Equation (44) is transformed in a set of real equations, by introducing the following definitions: where the phases are set to make the system autonomous as: Then, real and imaginary parts of the equations are collected, giving rise to the following set of six firstorder differential equations in the variables x k , y k : Equilibrium points of Eq. (47) are sought, and their stability is analyzed by evaluating the eigenvalues of the corresponding Jacobian matrix. The real amplitudes are then evaluated as: whereas the motion of the system is reconstituted with Eqs. (36) and (41).

Case 2: Ω ω 3
For the specified case, the differences with respect to case 1 are highlighted. In particular, the following scaling is adopted for the defined forcing terms: The external detuning here is defined so that: while the internal detunings are those given in Eq. (31). The perturbation equations are the same as in the previous case, namely Eqs. (32)-(34), where the forcing terms is now:

Linear-order problem
The solution of the linear-order problem is: where the components of the vectors Λ j are:

Quadratic-order problem
By substituting expressions (52) into Eq. (33), the second-order amplitude modulation equations reduce to: and the particular solution of the quadratic-order problem Eq. (33) is: Vectors w j appearing in Eq. (55) are not explicitly written in this specific case for the sake of brevity.

Cubic-order problem
By substituting expressions (52) and (55) into Eq. (34), the amplitude modulation equations at the third order are derived, namely: Reconstitution is made as in Eq. (44) and then, making use of Eq. (45), with the following definition of phases: the set of real ordinary differential equations is obtained: The numerical values of the coefficients appearing in Eq. (58) are explicitly given in Appendix E for the numerical application defined in the next section, whereas their analytic expressions are omitted for the sake of brevity.
Still, equilibrium points and relevant stability conditions are evaluated, and the real amplitudes r k = x 2 k + y 2 k analyzed (k = 1, 2, 3).

Numerical results
The following parameters are assumed for the pipe under analysis: mean radius of the cross-section R = 0.1 m, thickness h = 4 mm, Young modulus E = 1.65 · 10 8 Pa, modal damping factor ζ = 1% for all the modes, while the length is varied in the range l ∈ [1, 3] m. The first three natural frequencies, evaluated from the eigenvalue problem Eqs.  Therefore, for both the aforementioned cases, it is chosen to have almost perfect internal resonance between the global mode (frequency ω 1 ) and the local one which is in 1:1 resonance with the external force, namely ω 2 for case 1 and ω 3 for case 2. Furthermore, the next closest frequency (ω 4 ) is much larger than ω 3 , this justifying the choice of the reduced basis of three modes in the Galerkin projection. Anyway, a convergence analysis of the reduced system in terms of number o (Fig. 5)f involved modes is performed as well. More specifically, numerical integration of the system projected on the first ten modes, i.e., adding further four global and three local modes, is carried out (details of the formulation are omitted).

Case 1: Ω ω 2
In this case, the external force produces a primary resonance on the second mode (the first local), the latter being in internal resonance with the other two modes (the first global and second local).  The frequency response curves, shown in Fig. 6, are directly expressed in terms of r k , k = 1, 2, 3, and are represented versus the external frequency, normalized with respect to the first modal frequency: Ω/ω 1 ; the stable branches of the solution are represented by the solid lines, while the unstable branches by dashed lines. Moreover, they are reproduced for different values of the internal detuning ρ 2 , in order to analyze in more detail the effect of the internal resonance: ω 2 = 2ω 1 (dark blue lines), ω 2 = 2.01ω 1 (light blue lines), ω 2 = 1.99ω 1 (cyan lines). For the three internal detuning values, the second amplitude r 2 always exhibits the typical behavior of a monomodal solution which becomes unstable in between two bifurcation points (Ω/ω 1 1.97, 2.02), where the parametric resonance with the first mode takes place (see Fig. 6b). Accordingly, the amplitudes r 1 , r 3 are zero outside the range in which the parametric resonance is activated, whereas the solution becomes tri-modal inside (see Fig. 6a, c). The effect of the internal detuning ρ 2 is to slightly distort the curves as well as slightly shift them toward lower values of Ω/ω 1 as ρ 2 is increased; however, it can be concluded that ρ 2 qualitatively leaves the phenomena essentially unchanged.
To validate the results derived via the perturbation solution, a comparison of the reconstituted solution in terms of peak values of the modal amplitudes q 1 , q 2 , q 3 is made with the outcomes of numerical integration of Eq. (25) carried out via a Runge-Kutta routine in Mathematica [27]. The comparison between approaches is conduced only in the case ω 2 = 2ω 1 , and it is illustrated in Fig. 7, where the blue lines indicate the perturbation solution, while the black dots denote the numerical results. It can be observed that the stable branches of the solution are very well captured for what concerns q 1 (see Fig. 7a) and q 2 (see Fig. 7b). However, the accuracy reached for q 3 , which by the way assumes much lower values than q 1 and q 2 , is slightly inferior (see Fig. 7c), perhaps due to the internal detuning ρ 3 which is quite large, although the behavior is still qualitatively quite well-captured.
The results are also compared in terms of time histories of modal coordinates, evaluated at Ω/ω 1 = 1.99. Those are reported in Fig. 8, where the blue lines denote the perturbations solution, while the numerical results are represented by the black lines. Specifically, as illustrated in Fig. 8a, q 1 is very well captured by the perturbation solution that completely overlaps the numerical solution. Similarly, q 2 is very well captured by the per- vs Ω/ω 1 ; c max(q 3 ) vs Ω/ω 1 . Solid line: perturbation method; dotted line: numerical integration turbation solution that overlaps the numerical solution except for a negligible difference (less than 2%) in correspondence of the peaks (see Fig. 8b). On the other hand, as for the frequency plots, the time history for q 3 highlights a slight loss of accuracy (see Fig. 8c). For that, the fitting may be improved by considering higherorder terms in the perturbation solution; however, as it is deduced by Fig. 9, the error in q 3 does not significantly affect the response in terms of global (v, ϑ) and local (a p , a w ) displacement variables of the beam-like structure. In particular, in Fig. 9, the peak response is evaluated at the beam tip s = l for the displacement v (see Fig. 9a) and the cross-section rotation ϑ (see Fig. 9b), and at s = l/4 for flattening a p (see Fig. 9c) and warping a p (see Fig. 9d) amplitudes. The numerical response of v and ϑ, strictly related to q 1 (see Eq. (24)), exhibits a very good agreement with the perturbation solution. Moreover, the response of the local variables a p (see Fig. 9c) and a w (see Fig. 9d) is strongly influ- Fig. 9 Frequency response curves for l = 1.246 m, ω 2 = 2ω 1 , ω 3 = 2.8ω 1 and Ω = ω 2 + 2 σ : a max(v(l)) vs Ω/ω 1 ; b max(ϑ(l)) vs Ω/ω 1 ; c max(a p (l/4)) vs Ω/ω 1 ; d max(a w (l/4)) vs Ω/ω 1 . Solid line: perturbation method; dotted line: numerical integration of the 3-d.o.f. system; crosses: numerical integration of the 10-d.o.f. system enced by q 2 , and the agreement between the perturbation solution and the numerical result is very good as well. This reveals that the contribution of the third mode is small, though not negligible since it turns out to have a significant role in the determination of the position of the bifurcation points. Furthermore, the grey crosses, indicating the outcomes given by integration of the ten d.o.f. system, provided for convergence analysis, show good qualitative agreement, with a small quantitative deviation in terms of amplitudes of limit cycles close to the right bifurcation point. The latter aspect confirms the validity of the three-mode reduction proposed here.
However, as a major result, Fig. 9 proves the energy exchange from the directly excited local to the global response of the pipe, due to the internal resonance.

Case 2:
In this case, the external force produces a primary resonance on the third mode (the second local), the latter being in internal resonance with the other two modes (the first global and first local).
The frequency response curves are shown in Fig. 10, directly expressed in terms of r k , k = 1, 2, 3. From this figure, it is clear that only r 3 is involved in the response, being the contribution of r 1 , r 2 vanishing. Therefore, in the considered case, even though internal resonances are present, the following two circumstances concurrently happen: (1) a de facto nonlinear orthogonality among modes occurs [28] in the selected range of frequencies, which induces the coefficients responsible for the modal interaction in the amplitude modulation equations (AMEs), even not zero, not to provoke coupling; (2) the nonlinear terms are not able to (independently) trigger the 1:3 and 2:3 subharmonic resonances on modes 1 and 2, respectively. As a consequence, the first global and local modes, besides the contribute induced by the external force, only passively participate to the motion. This occurrence is confirmed by the frequency response curves in terms of reconstructed modal coordinates q 1 , q 2 , q 3 , shown in Fig. 11. There, the solution obtained by the perturbation method, shown in blue solid lines, is superimposed to numerical solutions (dots) obtained by numerical integration of the Galerkin equations. On the one hand, as expected, the response has prevailing component on the amplitude q 3 , that is, directly activated by the external excitation. It exhibits a hardening behavior characterized by the presence of multi-valued solutions and unstable branch (see Fig. 11c). On the other hand, the response of q 1 , q 2 (see Fig. 11a, b respectively) is mainly governed by the non-resonant terms ensuing at the first and second order (see Eqs. (52) and (55)). It is observed that though q 1 , q 3 are very well captured, the perturbation solution slightly loses accuracy in terms of q 2 around Ω/ω 1 3.04 (see Fig. 11b).

Fig. 11
Frequency response curves for l = 1.344 m, ω 2 = 2.3ω 1 , ω 3 = 3ω 1 and Ω = ω 3 + 2 σ : a max(q 1 ) vs Ω/ω 1 ; b max(q 2 ) vs Ω/ω 1 ; c max(q 3 ) vs Ω/ω 1 . Solid line: perturbation method; dotted line: numerical integration As done for case 1, the response is also compared in terms of time histories of the modal coordinates evaluated at Ω/ω 1 = 3.04, and the result is illustrated in Fig. 12. As expected, the time histories of q 1 , q 3 predicted by the perturbation solution completely overlap the numerical curves (see Fig. 12a, c), while q 2 is affected by a slight loss of accuracy in correspondence of the peaks (see Fig. 12b).
Finally, the global and local variables response curves are reconstituted and compared to the numerical solution. As done previously, the maximum value of the response of v and ϑ are evaluated at the beam tip s = l, and the curves are represented in Fig. 13a, b, respectively, whereas the response of a p and a w is evaluated at s = l/4, and the curves are shown in Fig. 13c, d, respectively. The response of the global variables is led by q 1 , whereas the local variables mainly follow q 3 , having q 2 a lower contribution in qualitatively determining the overall response, though it has a significant role in quantitative terms. It can be finally observed that no exchange of energy from the local to the global Fig. 13 Frequency response curves for l = 1.344 m, ω 2 = 2.3ω 1 , ω 3 = 3ω 1 and Ω = ω 3 + 2 σ : a max(v(l)) vs Ω/ω 1 ; b max(ϑ(l)) vs Ω/ω 1 ; c max(a p (l/4)) vs Ω/ω 1 ; d max(a w (l/4)) vs Ω/ω 1 . Solid line: perturbation method; dotted line: numerical integration of the 3-d.o.f. system; crosses: numerical integration of the 10-d.o.f. system modes of the pipe occurs in this case. This is also confirmed by the outcomes of the 10 d.o.f. system, indicated by grey crosses; quantitative differences with the 3 d.o.f. system occur only on the passive mode amplitudes, which actually give a very small contribution to the overall response of the pipe, still confirming the validity of the three-mode reduction.

Conclusions
The nonlinear dynamic response of a pipe under external harmonic load is addressed in the paper. The pipe is modeled as a beam-like structure, taking into account the change in shape of the cross sections by means of the introduction of specific (local) variables. Specifically, the change in shape provides a further contribution to the system, which is competing with those given by elastic and inertial terms, which typically interact in cantilevers. The load, acting on half cap of the pipe, has nonzero direct component in the equation ruling local variables, and it is assumed resonant with one of the local modes. Moreover, 1:2 and 1:3 internal resonances between global and local modes are considered as well.
After a Galerkin projection, the response of the pipe is addressed for two different load cases, respectively, i.e., external resonance with the first or second local mode, via perturbation methods. Specific scaling and implementations of the MSM are carried out for the two cases.
The obtained solutions, which are in general good agreement with numerical integration, show possible exchange of energy from the local to the global motions. On the one hand, if the external force is resonant with the first local mode, the nonlinear terms are able to trigger the internal resonances and induce tri-modal solutions. On the other hand, if the load is resonant with the second local mode, the internal resonances are not activated due to a de facto nonlinear orthogonality among modes, and the response mostly remains bounded in the local behavior.
Author contributions All authors contributed to the study conception and design. Material preparation, data collection, and analysis were performed by Arnaldo Casalotti and Daniele Zulli. The first draft of the manuscript was written by Arnaldo Casalotti and Daniele Zulli, and all authors commented on previous versions of the manuscript. All authors read and approved the final manuscript.
Funding Open access funding provided by Università degli Studi dell'aquila within the CRUI-CARE Agreement. The authors declare that no funds, grants, or other support were received during the preparation of this manuscript.

Data availability
The datasets generated during and/or analyzed during the current study are available from the corresponding author on request.

D Terms of the particular solution
The column vectors appearing in the particular solution at the second order (Eq. (41)) for the case 1 assume the following expressions: The numerical values of the coefficients appearing in the amplitude modulation equation for case 1 (Eq. (47)) are given as follows, with reference to the case study of Sect. 5.1: