Nonlinear dynamics of traveling beam with longitudinally varying axial tension and variable velocity under parametric and internal resonances

In this investigation, the nonlinear dynamics of an axially accelerating viscoelastic beam on a pully mounting system have been analyzed. The axial tension of the beam is modeled as a function of the traveling velocity, support stiffness parameter as well as spatial coordinate. Geometric cubic nonlinearity in the equation of motion is due to elongation in the neutral axis of the beam. The integro-partial differential equation of motion of the axially accelerating beam associated with the simply supported end conditions is solved analytically by adopting the direct perturbation method of multiple time scales. As a result, a set of complex variable modulation equations is generated, which governs the modulation of amplitude and phase. This set of modulated equations is numerically solved to explore the influence of the support stiffness parameter upon the stability and bifurcation of the beam which has not been addressed in the existing literature. Apart from this, the impact of fluctuating velocity component, viscoelastic coefficient, longitudinal stiffness parameter, internal and parametric frequency detuning parameters on the stability and bifurcation analysis is studied, revealing significant dynamic characteristics of the traveling system. The fourth-order Runge–Kutta method is applied is to find the dynamic solution of the system. The system displays stable periodic, quasi-periodic, and mixed-mode dynamic responses along with the unstable chaotic behavior for a specific set of system parameters. The results obtained through an analytical–numerical approach may help the design and operation of an axially moving beam.

the equation of motion of an axially moving system. Pakdemirli and Ozkaya [2] and Yang et al. [3] have examined the behavior of traveling continua that maintained constant velocity. This axial velocity is modeled as a harmonic function in the investigations carried out in [4][5][6][7][8][9][10][11][12][13][14][15][16] by considering different beam materials; damping terms; linear and nonlinear models of the beam. Oz [4], Pakdemirli and Oz [5], and Oz et al. [6] used elastic beam material with simply supported boundary conditions. Pakdemirli and Oz [5] introduced geometric nonlinearity into the system. The effect of damping on traveling beams is extensively investigated considering the linear beam model [6] and nonlinear beam model [7]. Marynowski et al. [8,9] studied the effect of internal damping focusing on different viscoelastic material models. They considered the Kelvin-Voigt model and Burgers models in [8] for a traveling web and concluded that for lower values of damping coefficient, both cases confer similar dynamic results, however for higher values of the same damping, the later model gives more appropriate results. Moreover, Marynowski [9] carried out similar studies on traveling viscoelastic beams considering Kelvin-Voigt and Maxwell models and obtained different numerical results at higher values of damping coefficient. Lv et al. [10] studied the traveling viscoelastic sandwich beam using Galerkin truncation followed by the method of multiple scales considering parametric excitation originating from variable speed. Saksa et al. [11] compared the dynamic characteristics of steel and paper adopting Poynting-Thomson model of axially accelerating viscoelastic beam. Yao and Yimin [12] carried out reliability and sensitivity analysis on traveling beam and concluded that the reliability of the moving system decreases when traveling velocity increases, the length of the beam increases for a particular speed, and it increases when the thickness of the beam increases and mass density decreases. Chakraborty and Mallik [13] made free and forced nonlinear vibration analysis of a traveling system using complex normal modes. Mao et al. [14] investigated Euler-Bernoulli viscoelastic moving beam excited parametrically. The multiple scale method (MMS) is used to reduce the higherorder equations into a set of first-order equations. It is observed from numerical solutions that the amplitudefrequency response plots tilt toward the left, due to the presence of quadratic nonlinearity in the system. Ding and Chen [15] examined the nonlinear dynamic characteristics of viscoelastic traveling beams employing the differential quadrature method (DQM). Ghayesh et al. [16] derived the nonlinear equation of motion in transverse and longitudinal directions adopting the force and moment balance approach and solved numerically by adopting a variable step-size modified Rosenbrock method. Wang et al. [17] focused on the axially accelerating hyper-elastic beams using MMS to solve the higherorder partial differential equation of motion analytically. Further, they applied Galerkin method for verification of numerical results obtained from the MMS approach. Padioussis and Moon [18] examined the dynamical behavior of a pipe carrying fluid with fixed-free support conditions. Whereas Czerwinski et al. [19] investigated pulsating fluid flows through a slightly bent pipe with clamped-clamped end conditions subjected to combination parametric resonance.
The above studies [2][3][4][5][6][7][8][9][10][11][12][13][14][15][16][17][18][19] have presented a large number of numerical examples exclusive of internal resonance. Several authors have paid attention to consider the internal resonance in the traveling systems with harmonically varying speed to get rich and interesting dynamics. Internal resonance carries theoretical importance since this phenomenon involves energy exchange between two different modes of oscillation in a continuous system. Few selected relevant studies [20][21][22][23][24][25][26][27][28][29] are presented where internal resonance is included in the investigation of different types of beam models using various approaches under different parametric resonant conditions. Chin et al. [20] analyzed the nonlinear response of moving beam with hinged-clamped end conditions subjected to different parametric resonant conditions. Huang et al. [21,22] investigated the nonlinear vibration of axially moving elastic beam excited with the transverse harmonic load. The dynamic behavior is evaluated using the incremental harmonic balanced (IHB) method [21]. Moreover, the same problem is extended in [22] by employing a new incremental harmonic balanced (IHB) method with two-time scales. In both the studies, the results obtained from the above-mentioned methods give an excellent match with numerically integrated solutions displayed in phase plane plot, time history, Poincare sections, Fourier spectra, and Lyapunov exponent. Ghayesh et al. [23] analyzed the nonlinear dynamic response of moving systems with time-varying axial velocity. The Galerkin technique is used to reduce the differential equation of motion into a set of ODEs. These ODEs are solved numerically using the Pseudoarclength method to plot the frequency response curve of the system. Mao et al. [24] examined the traveling beam under parametric excitations and 3:1 internal resonance. When the system is excited parametrically with the first mode, the zone of unstable trivial response in both modes are affected by internal resonance and saddle type bifurcation point comes into the picture; but when the parametric excitation of the second mode is considered, only the second mode exhibits the parametric responses. However, under combination parametric resonance, the system displays different results compared to the previous two cases.  studied the various aspects of traveling beams. They show the effect of viscoelastic behavior and internal resonance upon the standard linear model of the beam in [25]. For the first time, they used multi-scale method with the real modal function to determine the steady-state response of the system; considered the effect of rotary inertia on the free vibration of moving beams with two different kinds of support conditions in [26], and introduced the nonhomogeneous boundary conditions into the nonlinear system in [27]. Wang et al. [28] studied the viscoelastic beam moving under parametric excitation in presence of 2:1 internal resonance. The third-order perturbation solution is used to examine the effect of quadratic and cubic nonlinearities more precisely. The frequency-response and force-response plots are used to illustrate the stability of the system. Wang et al. [29] studied the nonlinear dynamic behavior of a system moving in a constant magnetic field experiencing 1:3 internal resonance.
The above-mentioned investigations focused on axially moving systems without internal resonance [2][3][4][5][6][7][8][9][10][11][12][13][14][15][16][17][18][19] and with internal resonance [20][21][22][23][24][25][26][27][28][29] with one thing common, i.e., the axial tension which remains constant throughout the length of the beam. If the axial tensions are of equal magnitude and opposite in direction at both ends of the moving system, then the beam moves with zero acceleration. But in reality, achieving a constant axial tension is quite difficult due to the presence of manufacturing imperfections, wear, and tears in the moving systems. There are few works [30][31][32][33][34][35][36][37] in which variable axial tension is assumed in axial traveling continua. Mote [30] examined the band saw vibration by considering the longitudinal tension of the band saw depending upon band velocity and pully mounting system. It was found that the natural frequency of the band drops off uninterruptedly when the band velocity increases to a higher value. Mockensturm et al. [31] studied the case of the viscoelastic translating string focusing on the dynamic tension that develops in the string when the system is stretched physically. Marynowski et al. [32] analyzed the dynamic characteristic of a viscoelastic moving beam with the Zener rheological model subjected to timevarying axial tension. Ma et al. [33] introduced the nonhomogeneous parabolic tension profile at either end of the traveling web and observed the decrease in natural frequency; critical velocity; and the change in mode shape of the moving web. Wang et al. [34] investigated the hyper-elastic beam subjected to timevarying axial load along with 2:1 internal resonance. In this work, the eigenfunction expansion method is employed analytically to get a series of nonlinear ODEs which are solved by the harmonic balance method followed by the pseudo-arc-length numerical approach to plot the amplitude-frequency response diagrams. Guo et al. [35] studied the dynamical behavior of composite pipe conveying fluid excited with harmonical axial tension. The numerical examples show the influence of axial tension and fiber orientation angle on the natural frequencies of the system. Lenci et al. [36] investigated the linear composite beam having two elastic layers of the Euler-Bernoulli type, and the influence of various boundary conditions in detail. Furthermore, Lenci and Clementi [37] examined the free vibrational analysis of a moving beam subjected to arbitrarily distributed axial load, which when attained a fixed value led to 1:3 internal resonance occurs between the first two modes. Furthermore, Lenci et al. [38] compared the nonlinear characteristics of the Timoshenko beam by incorporating mechanical and geometrical curvatures. A relatively more appropriate approach is to assume the velocity as well as tension to be variable. Few works [39][40][41][42][43][44][45][46][47] are available to fill up this lack of research gap for a nonlinear traveling system. Chen and tang [39][40][41] investigated the vibrational analysis of axially accelerating viscoelastic beam with nonhomogeneous boundary conditions. The said system is excited parametrically due to longitudinally varying tension and harmonically varying velocity [39]. Further, Tang and co-workers carried out similar studies on moving beams by introducing 3:1 internal resonance [40] and considering the velocity and longitudinal tension as a function of time and interdependence [41]. Zhang et al. [42] analyzed the effect of two frequency excitations on the stability boundary of the traveling beam. Lv et al. [43] examined axially moving viscoelastic sandwich beam where the axial tension and velocity are considered as harmonic functions of time considering the Euler beam model [43] and Timoshenko beam model [44]. Yan et al. [45] studied the effect of longitudinal varying tension and time-varying velocity on functionally graded beam material. Tang et al. [46,47] examined the nonlinear dynamics of the traveling beams by adopting the method of multiple time scale (MMS) followed by the differential quadrature method (DQM). The influence of velocity-dependent tension and tension-dependent velocity are investigated in [46] while the time and space-dependent axial tension and tension-dependent axial speed are analyzed in [47]. Both the studies concluded that with increasing damping parameters the instability region appearing between the different bifurcation points decreased.
From the literature survey, it has been observed that there is no detailed study on stability and bifurcation analysis considering the effect of support stiffness of a nonlinear traveling beam under the effect of simultaneous parametric and internal resonant conditions. This study investigates the nonlinear planer vibration of the viscoelastic beam moving over the elastic support with simply supported boundary conditions. The frequency of parametric excitation is equal to two times of first natural frequency of the system in presence of 3:1 internal resonance. The longitudinal tension is considered as a function of velocity, support stiffness parameter, and it varies over the full span of the beam. At the same time, the axial speed is assumed to be harmonically changing over a mean value. To solve the cubic nonlinear integro-partial differential equation of motion, the perturbation method, i.e., the Method of Multiple time Scales (MMS) is adopted and a set of complex variable modulation equations results. These modulation equations are analyzed numerically for stability and bifurcations of the system with different values of system parameters like pulley stiffness parameter, fluctuating velocity component, viscoelastic coefficient, longitudinal stiffness parameter, internal and parametric detuning parameters. Also, dynamic responses of the system are studied to illustrate a wide variety of dynamic behavior like periodic, quasiperiodic, mix-mode, and chaotic responses.

Formulation of the problem
A uniform horizontal beam of simply supported end conditions is considered as shown in Fig. 1. The assumptions made are (1) velocity is varying harmonically over a mean value and axial tension varies longitudinally over the full length of the beam (2) the beam is experiencing planner motion (3) the transverse plane of the beam remains plane during the motion and the beam behaves like Euler-Bernoulli beam in transverse vibration (4) geometric cubic nonlinearity arises due to stretching of the neutral axis of the beam (5) the lateral motion of the beam is not considered.
The equation of motion may be reproduced by incorporating the cubic nonlinearity due to mid-plane stretching [1,22], material damping [16,26,27], and external viscous damping [13,25] in the form With boundary conditions In the above Eqs. (1-2), the star mark denotes that the above parameters are dimensional ones, u denotes the deflection in the transverse direction of the beam. V and P denote the axial velocity and the axial tension of the beam, respectively. q, m , A , L , C , E denote the mass density, mass per unit length, cross-sectional area, length, the external damping factor, and the coefficient of Young's modulus of the beam, respectively.
The variable axial tension of the traveling beam may be expressed as [39,40,42] Equation (3) can be interpreted with practical as well as theoretical significance. In most of the studies , it has been noticed that magnitude of the longitudinal tension of the moving systems was considered to be independent of velocity and spatial coordinate. This signifies that the axial tension of the moving system is having equal magnitude and opposite direction at either end. This approach does not consider the fact that the traveling system moves with nonzero acceleration. Because the moving system has a time-dependent speed component (i.e., V(t) = V 0-? eV 1 sinXt), it develops nonzero additional inertial tension force as a consequence of Newton's 2nd law of motion. Many authors have not considered this aspect only to reduce mathematical complexity which is easy to analyze. Therefore, an appropriate model for the axial tension of the moving beam is expressed in Eq. (3). The first term is the mean axial tension; the second term signifies the dependence of tension on axial speed and also on the stiffness of the pulley mounting system g ¼ 1=ð1 þ k r =ð2EA=LÞÞ ½ where k r is the relative support stiffness parameter. Support stiffness parameter g varies between 0 and 1 [30]. g ¼ 1 when the system drives are mounted over the pulleys having no rigidity and g ¼ 0 when the pulleys have infinite rigidity. The third term ðx Ã À LÞqA oV Ã =ot Ã ð Þ stands for variation of axial tension along the longitudinal direction of the beam because of nonzero acceleration.
To convert the dimensional parameter to nondimensional one, the following scheme is used Applying axial tension Eq. (3) and nondimensional scheme Eq. (4) in the dimensional equation of motion Eq. (1) and boundary condition Eq. (2) and also considering the transformation u Ã ¼ ffiffi e p u Ã [5], weekly nonlinear and nondimensional equation of motion with the associated boundary condition may be obtained as where v f and v l represents the non-dimensional flexural stiffness and longitudinal stiffness of the beam respectively, l is the nondimensional external damping factor, a is the dimensionless material damping parameter, and e is the bookkeeping parameter e\1 ð Þ. The subscripts t represents the differentiation with respect to time and the subscript x denotes the derivative with respect to spatial coordinate.
The traveling velocity of the beam is assumed as a harmonic function.
where V 0 is mean axial velocity, eV 1 is the amplitude of fluctuating velocity, and X is the frequency of fluctuation. This assumption can have a physical interpretation. For example, when the moving beam (for example a belt) has planner motion over a pair of pulleys, the torsional oscillation of the pulley can lead to set a small fluctuation in the axial speed. The variation of the axial speed may also be generated from other sources, for example, due to eccentricity, wear and tear of the pulley. The speed fluctuation induced due to these factors may be periodic but having a major Fourier harmonic component eV 1 sin X t, which excites the system parametrically and develops various parametric resonant states in the system. Substituting Eq. (7) in Eq. (5), we get Fig. 1 The schematic of the traveling beam with simply supported end conditions With the boundary conditions Equation (8) contains a cubic geometric nonlinearity term, which does not lead to closed form of solution. Hence to get an approximate solution a perturbation method is used.

Method of analysis
To get an approximate solution for nonlinear system (Eq. 8), a perturbation method, i.e., the method of multiple time scale (MMS) is adopted [4,5,20]. In the direct perturbation method, one can consider the expansion in the form of where u 0 and u 1 are the zeroth-order and first-order terms, the time scale used here is T n ¼ e n t where n ¼ 0; 1; 2; . . . for n = 0, T 0 = t represents the fast time scale and for n = 1, T 1 ¼ eT 0 represents the slow time scale. The first-and second-order derivatives of the time scale in terms of T 0 and T 1 are By using Eqs. (10) and (11) in the partial differential Eqs. (8)(9) and separating the coefficients of like powers of e, we get The solution of Eq. (12) may be expressed in this manner where x n is the natural frequencies of the continuous system and / n is the complex mode shapes, which may be presented as below [4] where b in are the Eigenvalues of the system which satisfies the dispersive relation Eq. (16) and end conditions Eq.
In this study, the principal parametric resonance of first mode X % 2x 1 ð Þis considered, i.e., the excitation frequency is equal to two times of the first natural frequency. With increasing mean velocity of the moving system, the centrifugal force reduces the flexural stiffness as a result the linear natural frequencies of the traveling beam decreases. For the specific combination of system parameters, the first and second modes of natural frequency have a commensurate relationship. This leads to internal resonance between the first two modes. From Table 1 the bold number row provides the information about the 3:1 internal resonance which takes place between the first two modes Þfor a specific value of the mean velocity V 0 ¼ 0.5145. All higher modes except the first and second modes decay with time due to the presence of damping terms in the equation of motion, thus first two modes will give the long-term response to the system so that we can replace Eq. (14) in the form The resonant frequency relations are Here r 1 and r 2 denote the frequency detuning parameters to determine the system behavior in the vicinity of the resonant frequencies of the traveling beam. By substituting Eqs. (18)(19) into Eq. (13), one can get where 1 1; 1 2; . . .1 9 are presented in the Appendix. The NST stands for the nonsecular terms which don't give long-term responses of the system. The righthand side of Eq. (20) corresponding to the nonhomogeneous part has a solution only if a solvability condition is satisfied [48], whereas the left-hand side of the equation corresponds to the homogeneous part with the associated boundary condition having a nontrivial solution. This requires the nonhomogeneous part of Eq. (20) to be orthogonal to every solution of the adjoint homogeneous problem, which leads to two complex variable modulation equations for phase and amplitude as given below in line with [20] The prime represents the differentiation with respect to time and S i ; g i ; C i ; k i ; e i are presented in Appendix. Over bar signifies the complex conjugate.

Stability and bifurcations analysis
From the complex variable modulated Eqs. (21)(22), the equilibrium solutions are developed and their stability and bifurcation analysis is executed for the case of X % 2x 1 . For this, Cartesian form of transformation is used for the complex amplitude A n as Substituting Eq. (23) into Eqs. (21) and (22), separating the real and imaginary parts after due simplification through standard mathematics, one may arrive at the normalized reduced equation or, the Cartesian form of modulation equations as given below Where k 0 1 ¼ 0:5r 2 and k 0 2 ¼ 1: To evaluate the stability of the system, the above equations are perturbed to obtain where T stands for transpose, the stability and bifurcation of the system are determined by the eigenvalues of the Jacobian matrix J c ½ . By setting up p 0 (24)(25)(26)(27) the steady-state response of the system is obtained.

Results and discussions
The steady-state response of the system is determined by solving the normalized reduced Eqs. (24-27) after setting continuation algorithm based on Pseudo-arc-length continuation technique [49]. At each point of the equilibrium solution, the eigenvalue of the Jacobian matrix is evaluated to determine the stability and bifurcation of the nonlinear system. The frequency response curves are symmetric about the r 2 axis, so the positive side of the response curves is shown on the frequency response diagram.

Frequency response analysis
Frequency response curves of the first and second modes are plotted between the variation of parametric frequency detuning parameter r 2 and amplitude of oscillation adopting continuation algorithm [49] which describes the stability and bifurcation of the system by considering system parameters V 1 ¼ 5; r 1 ¼ 68:2; k ¼ 0:99; a ¼ 0; l ¼ 0:1; v l ¼ 40; and v f ¼ 0:2; as shown in Fig. 2. All the bifurcation branches of frequency response curve are tilted toward the right and this signifies that the system exhibits the hardening type of nonlinearity. The continuous line represents the stable equilibrium solutions, the bold line indicates the unstable foci, and the dotted line indicates the saddle type instability. When the value The amplitude of this branch (P 1 -P 3 ) rises initially and then falls in the first mode (Fig. 2a), whereas in the second mode (Fig. 2b)   . From frequency response plot, it has been observed that, the amplitude of isolated curve (E-F-G) reduces in the first mode, whereas its opposite trend is obtained in second mode. Furthermore, decreasing r 2 along the same isolated curve the equilibrium solution meets two more Hopf bifurcation points, i.e., H 6 ðr 2 ¼ 127:2527Þ and H 7 ðr 2 ¼ 85:1527Þ the region between these two points is unstable. Moreover, decreasing the value of r 2. Along the trivial branch (B-A), the trivial equilibrium solution passes over the Hopf bifurcation type unstable zone, i.e., H 12 À H 11 j j¼ 110:4000 and a reverse pitchfork bifurcation point P = 2 ðr 2 ¼ 45:2644Þ.

Effect of material damping parameter (a)
This section studies the impact of material damping parameter on stability and bifurcation analysis curves by employing the different value of material damping parameter, i.e., a ¼ 0:0001, a ¼ 0:0005 and a ¼ 0:001 the results are presented in Figs. 3a, b, 4a, b, and 5a, b respectively by keeping other system parameters as v f ¼ 0:2, v l ¼ 40, k ¼ 0:99, r 1 ¼ 68:2, l ¼ 0:1, V 1 ¼ 5. By comparing Figs. 2, 3, 4, 5 some qualitative and quantitative changes have been observed and summarized in Tables 2, 3. When the value of material damping parameter increases the nontrivial isolated curve (E-F-G) vanishes in case of a ¼ 0:001 Fig. 5a, b. Apart from this, the number of pitchfork bifurcation pointes are reduced to two viz. P 1 (supercritical pitchfork bifurcation point) and P 2 (reverse pitchfork bifurcation point) corresponding to the case of a ¼ 0:0005 (Fig. 4) and a ¼ 0:001 (Fig. 5), while it is unaltered associated with case of a ¼ 0:0001 (Fig. 3). Increasing r 2 along the nontrivial bifurcation branch ðP 1 À P 3 Þ, the number of bifurcation points are reduces and settle down to two, i.e., ðH 1 and SN 1 Þ in case of a ¼ 0:0001 ( Fig. 3) whereas, it is remain same in case of a ¼ 0:0005 (Fig. 4) comparing with the previous case a ¼ 0 (Fig. 2). Moreover, this nontrivial branch follows another trend Þwith the presence of bifurcation points, viz, ðH 1 ; H 2 ; H 3 ; SN 1 ; H 4 ; H 5 ; SN 2 ; H 6 ; and H 7 Þ related to the case of a ¼ 0:001 (Fig. 5a, b). Two extra saddle-type bifurcation points SN 4 ð Þ and SN 5 ð Þincorporates with the nontrivial branch ðP 2 À P 3 Þ in case of a ¼ 0:0005 (Fig. 4) and two additional Hopf bifurcation points viz. H 13 and H 14 are appeared on trivial branch (A-B) corresponding to the case of a ¼ 0:0001 (Fig. 3). If parametric detuning parameter r 2 decreases from a higher value through the isolated curve (E-F-G), the stable region widens significantly from H 4 À H 5 j j a¼0 ¼ 35:5201 to H 4 À H 5 j j a¼0:0001 ¼ 216:7124 with rising the influence of a ¼ 0 (Fig. 2) to 0.0001 (Fig. 3), furthermore, increasing a to 0.0005 (Fig. 4) the instability region reduces significantly and appears between the said points H 4 À H 5 j j a¼0:0005 ¼ 45:5193. From Tables 2, 3 it is found that the saddletype unstable portion between ðP 1 À P 2 Þ almost same for all cases of material damping parameter along the trivial branch (A-B) Fig. 7a, b the pattern of the said nontrivial branch has modified to P 1 À C ð Þ     branch with the presence of three bifurcation points like ðH 1 ; SN 1 ; and SN 2 Þ. With increasing r 2 along the trivial branch (A-B), the trivial equilibrium solution meets a reverse pitchfork bifurcation point P 2 where the nontrivial bifurcation branch P 2 À P 3 ð Þoriginates in Fig. 6a, b and P 2 À E ð Þin Fig. 7a, b. Corresponding to the P 2 À P 3 ð Þbranch two additional saddle-node bifurcation points are appearing at P 6 ðSN 4 Þðr 2 ¼ 56:2927Þ and P 7 ðSN 5 Þðr 2 ¼ 54:4427Þ shown in (Fig. 6), which is missing in the earlier case of k = 0.99 (Fig. 2) also in case of k = 0.5 (Fig. 7). Furthermore, increasing r 2 along the same trivial line ðA À BÞ, the equilibrium solution approaches to Hopf bifurcation point H 8 and supercritical pitchfork bifurcation point P  Tables 4, 5. Decreasing r 2 from a higher value, the nontrivial equilibrium solution either follows to the isolated curve E À F À G ð Þin case of k = 0.99 (Fig. 2) and k = 0.7 (Fig. 6), whereas the said isolated curve is disappears in case of k = 0.5 (Fig. 7) and it follows to the nontrivial branch C À P 1 ð Þ or follows to trivial line B À A ð Þ. The maximum portion of this E À F À G ð Þ curve and C À P 1 ð Þbranch goes through unstable zone up to the saddle-node bifurcation point P 4

Effect of external damping parameter (l)
A through study on the impact of variation of external damping parameter on stability and bifurcation analysis is carried out by keeping other system parameters same as before, i.e., V 1 ¼ 5; r 1 ¼ 68:2; k ¼ 0:99; a ¼ 0; v f ¼ 0:2; and v l ¼ 40. Different viscous medium is considered for the traveling beam like, l ¼ 0:1 (Fig. 2a, b), l ¼ 0:5 shown in (Fig. 8a,  b), and l ¼ 0:7 presented in (Fig. 9a, b) and the respective values of r 2 for all bifurcation points are abstracted in Tables 6, 7. From the above figures and tables some qualitative and quantitative similarity dissimilarities have been observed. The qualitative similarity comprises with the nontrivial isolated curve E À F À G ð Þ , which is present in all three cases of external damping parameter. The quantitative similarity contains the frequency response plots have equal number of pitchfork bifurcation points, i.e., ðP 1 ; P 2 ; P = 1 ; and P = 2 Þ corresponding to the cases of l ¼ 0:1 and 0:5, but it reduces to two, viz. ðP 1 and P 2 Þ in case of l ¼ 0:7 (Fig. 9a, b) which incorporates with the quantitative dissimilarity. Besides this, two additional saddle node bifurcation points like P 6 ðSN 4 Þ and P 7 ðSN 5 Þ are appears on the nontrivial branch ðP 2 À P 3 Þ associated with the case of l ¼ 0:7 (Fig. 9), which was absent in last two cases of Þ were present in earlier case l ¼ 0:1 (Fig. 2), are disappears from the present cases of l ¼ 0:5 and 0:7. Increasing parametric detuning parameter, along the nontrivial two-mode solution branch P 1 À P 3 ð Þ the unstable region reduces from    Fig. 11 Amplitude-plots a first mode, b second mode for r 2 ¼ 50;r 1 ¼ 68:2; k ¼ 0:99;a ¼ 0;l ¼ 0:1;v f ¼ 0:2; and v l ¼ 40 r 2 from a higher value it may approaches to isolated nontrivial curve E À F À G ð Þ , the Hopf bifurcation type instability portion decreases significantly from jH 4 À H 5 j l¼0:1 ¼ 35:5475 (stable region) to jH 4 À H 5 j l¼0:5 ¼ 66:3370 and jH 4 À H 5 j l¼0:7 ¼ 40:8910. Including this, the Hopf-related unstability region also compressed between jH 12 À H 11 j with increasing the l value like, from jH 12 À H 11 j l¼0:1 ¼ 110:4000 to jH 12 À H 11 j l¼0:5 ¼ 27:6774 and jH 12 À H 11 j l¼0:7 ¼ 21:8137.

Amplitude response analysis
The amplitude response plots shown in Fig. 10a, b describe the response of the system against the harmonically varying velocity component V 1 ð Þ considering the system parameters r 1 ¼ 68:2, r 2 ¼ 100, ð Þ , where the system response either jumps to the nontrivial branch P 1 À G ð Þor to the isolated nontrivial curve D À C À E ð Þas per the region of attraction and initial condition. Increasing V 1 through P 1 À G ð Þ branch it approaches two saddle node bifurcation points and one Hopf bifurcation point SN 1 ; H 2 ; and SN 2 respectively, the equilibrium solutions are stable between SN 1 V 1 ¼ 0:0369 ð Þ and H 2 V 1 ¼ 0:6400 ð Þ whereas it is unstable between H 2 V 1 ¼ 0:6400 ð Þ and SN 2 V 1 ¼ 5:0924 ð Þ . Furthermore, increasing V 1 on same branch it encounters one saddle-type and Hopf-Related unstable points Decreasing the amplitude of axial velocity V 1 along the isolated curve D À C À E ð Þ , it is observed that a saddle node bifurcation point SN 5 V 1 ¼ 15:6775 ð Þ occurs at C and the solution becomes unstable throughout the curve D À C ð Þ, amplitude of this branch decreases continuously in first mode (Fig. 10a) but it decreases first then increases in second mode (Fig. 10b). This branch D À C ð Þ also attached with another saddle-type unstable branch C À E ð Þ.

5.6
Effect of parametric frequency detuning parameter (r 2 ) Figure 11a, b and Table 8 show the impact of parametric frequency detuning parameter ðr 2 Þ on the amplitude response curves taking r 2 ¼ 50 and keeping other parameters same as earlier case, i.e., r 1 ¼ 68:2;k ¼ 0:99;v l ¼ 40;v f ¼ 0:2;a ¼ 0; and l ¼ 0:1. The significant change observed in the present case is that the nontrivial curve P 1 À G ð Þ which appeared for r 2 ¼ 100 (Fig. 10a, b) disappears for r 2 ¼ 50 (Fig. 11a, b). Apart from it, the properties of the isolated nontrivial curve D À C À E ð Þhas been modified with the inclusion of five bifurcation points ðH 4 ; SN 5 ; SN 6 ; H 5 ; and SN 7 Þ corresponding to Fig. 11a, b as compare to the appearance of a single saddle-node bifurcation point ðSN 5 Þ on the same curve corresponding to the case of r 2 ¼ 100 (Fig. 10a, b), which are summarized in Table 8. The range of instability due to Hopf bifurcation along the trivial equilibrium solution A À B ð Þ decreases from H 1 À P 1 j j r 2 ¼100 ¼ 19.9712 to H 1 À P 1 j j r 2 ¼50 ¼ 9.3260.

Effect of material damping parameter (a):
This paragraph describes the influence of material damping parameter a upon amplitude response plots by keeping all other system parameters same as first case r 2 ¼ 100,  Table 9. Observing Figs. 10, 12, 13 and Table 9 it is found that, the isolated nontrivial curve D À C À E ð Þ vanishes from the amplitude response plot in case of a ¼ 0:001 (Fig. 12a, b) whereas the property of this isolated curve is modified with the absence of SN 5 corresponding to the case of a ¼ 0:003 (Fig. 13a, b). Besides this, the properties of nontrivial branch Þhas changed with the disappearance of Hopf bifurcation point H 2 associated with the case of a ¼ 0:001 (Fig. 12a, b) whereas three additional bifurcation points like ðH 3 ; H 4 ; H 5 ; H 6 ; and H 7 Þ are appears on the same branch in case of a ¼ 0:003 (Fig. 13a, b). The amplitude of this this nontrivial branch P 1 À G ð Þ corresponding to a ¼ 0:003 increases monotonically in first mode but it decreases first then increases continuously in second mode. Comparing the trivial solution branch AÀB ð Þ of all cases of a, the Hopf bifurcation type of unstable zone reduces from H 1 À P 1 j j a¼0 ¼ 19:9712 to H 1 À P 1 j j a¼0:001 ¼ 18:5161 corresponding to Fig. 12a, b and to H 1 À P 1 j j a¼0:003 ¼ 17:2781 associated with Fig. 13a, b.

Effect of support stiffness parameter (k)
This section reveals the effect of support stiffness parameter (k) on amplitude response curves freezing all other system parameters r 2 ¼ 100, r 1 ¼ 68:2, a ¼ 0, l ¼ 0:1, v f ¼ 0:2, and v l ¼ 40. Three values of rigidity parameter k ¼ 0:5; 0:7 and 0:99 ð Þ are considered and resulting bifurcation points are displayed in Table 10. Unlike the impact of variation of k on the (a) (b) Fig. 12 Amplitude-plots a first mode, b second mode for r 2 ¼ 100; r 1 ¼ 68:2; k ¼ 0:99; a ¼ 0:001; l ¼ 0:1; v f ¼ 0:2; and v l ¼ 40 Table 9 Different bifurcation points on amplitude plot for r 1 ¼ 68:2; r 2 ¼ 100;k ¼ 0:99; l ¼ 0:1;v f ¼ 0:2; and v l ¼ 40     Tables 4, 5, the amplitude response curves exhibit more similarity so far as the nature of equilibrium solution curves and their corresponding bifurcation points are concerned. The amplitude response curves corresponding to k = 0.5 is shown in Fig. 14a, b and that of for k = 0.7 is dropped to avoid redundancy. It may be noted that the Hopf bifurcation point Þwhich is present in the case of k = 0.99, disappears in the subsequent cases of k = 0.7 and 0.5. If we observe the range of Hopf -related instability in the trivial and nontrivial solution curves, an increasing trend is noticed with the decrease in the value of k. Using the data displayed in Table 10, it may be shown that, SN 2 À SN 1 j j k¼0:7 ¼ 4:8923 and SN 2 À SN 1 j j k¼0:5 ¼ 4:9915. Apart from this, the   The influence of the external damping parameter (l) on amplitude response analysis keeping all other system parameters intact, i.e., r 2 ¼ 100, r 1 ¼ 68:2, k ¼ 0:99, a ¼ 0, v f ¼ 0:2, and v l ¼ 40 has been studied. Three different denser mediums have been assigned for the moving beam with the nondimensional viscous coefficients, l ¼ 0:1; 0:5 and 0:7 the amplitude response plots are displayed in Figs. 10a, b,  15a, b and 16a, b respectively. The bifurcation points appearing in the steady-state response of all three cases are presented in Table 11. If one compares the figures corresponding to above cases, then qualitative similarities and significant quantitative similarities as well as dissimilarities are surfaced. Qualitative similarity shows that the system displays similar trivial, nontrivial and nontrivial-isolated curves along with same kinds (Hopf, saddle-type and super-critical pitchfork) bifurcation points. Quantitative similarities include the maximum amplitude of both the modes remain same for all three cases, i.e., l ¼ 0:1 (Fig. 10a, b, l ¼ 0:5 (Fig. 15a, b), and l ¼ 0:7. Besides this, the system is having single super-critical pitchfork bifurcation point (P 1 ) in all three cases of l. The moving system experiences noteworthy quantitative dissimilarities with the variation of l: Three additional Hopf bifurcation points H 4 ; H 5 and H 6 ð Þappear with the present case (l ¼ 0:5 and 0:7) when compared with case l ¼ 0:1: Apart from it, the range of instability between the saddle-node bifurcation points ðSN 1 and SN 2 Þ reduces with increasing the viscous damping parameter. The numerical values are enumerated from Table 11, i.e., SN 2 À H 2 j j l¼0:1 ¼ 4:4524, H 4 À H 2 j j l¼0:5 ¼ 1:0877 and H 4 À H 2 j j l¼0:7 ¼ 0:7544. The similar trend is observed for the range of Hopfrelated instability between H 3 and SN 3 along the nontrivial curve P 1 À G ð Þ, i.e., H 3 À SN 3 j j l¼0:1 ¼ 10:5325, H 3 À SN 3 j j l¼0:5 ¼ 9:5000 and H 3 À SN 3 j j l¼0:7 ¼ 9:0000. Furthermore, if one observes the nontrivial isolated curve D À C À E ð Þ , the stable equilibrium solution exists between the Hopf bifurcation points H 5 and H 6 in case of l = 0.5 and 0.7 (Figs. 15, 16). This stability region shrinks with the decrease in l i.e., H 6 À H 5 j j l¼0:7 ¼ 4:3000; H 6 À H 5 j j l¼0:5 ¼ 2:1413 and finally vanishes for l ¼ 0:1. Along the trivial equilibrium solution curve, the position of the supercritical bifurcation point(P 1 ) is insignificantly varying in all three values of l as seen from Table 11. However, the Hopf bifurcation point H 1 shifts toward right with the increase in l value variety of dynamic solutions, some interesting results are presented. The dynamic solution of the nonlinear system is highly initial point sensitive. Figure 17a-e illustrates the system behavior corresponding to the initial point at r 2 ¼ 21:7034 nearer to the Hopf bifurcation point H 1 along the nontrivial branch ðP 1 À P 3 Þ corresponding to the frequency response plot (Fig. 3)for system parameters V 1 ¼ 5, r 1 ¼ 68:2, k ¼ 0:99, a ¼ 0:0001, l ¼ 0:1, v f ¼ 0:2 and v l ¼ 40.
The system displays mixed mode behaviour, i.e., periodic in the first mode and period-2 in the second mode which is illustrated in form of phase portrait (Fig. 17a, b), and time history (Fig. 17c, d).Taking another initial point r 2 ¼ 33:5034 along the same branch, the systems display quasi-periodic behavior in the first mode and quasiperiodic with the chaotic tendency in the second mode as described in phase portrait (Fig. 18a, b) and Poincare maps (Fig. 18c, d).
At another higher value of r 2 ¼ 86:1034 nearer to the saddle-type unstable point P 3 ðSN 1 Þ on the same curve ðP 1 À P 3 Þ; it has been observed that both the modes display chaotic behavior as presented in form of phase portrait (Fig. 19a, b), time history (Fig. 19c, d), Poincare map (Fig. 19e), and FFT power spectra (Fig. 19f).
When the effect of external damping parameter rises to l ¼ 0:35 keeping other system parameters intact, i.e., V 1 ¼ 5;r 1 ¼ 68:2;k ¼ 0:99;a ¼ 0;v f ¼ 0:2; and v l ¼ 40, the nonlinear system experiences mix mode response with period-2 and periode-3 in first and second modes respectively at r 2 ¼ 35:6799 as depicted from phase portrait and time plots (Fig. 23a-d). At another point r 2 ¼ 40:5799 the equilibrium solution exhibits period doubling phenomena in both the modes where repetition occurs in phase portrait (Fig. 24a, b) and time plots (Fig. 24c,  d). Further increasing detuning parameter along the same curve quasiperiodic behavior is noticed at r 2 ¼ 51:9201 as demonstrated in Poincare maps (Fig. 24e,  f). For another higher value of r 2 ¼ 54:4201 the chaotic response is detected through phase portrait, time plots, Poincare maps and FFT spectra as shown in Fig. 25a-f respectively. The transformation is detected in the system behavior from stable-periodic to chaotic through quasi-periodic route, for different combination of system parameters. The above-mentioned transition happens in the proximity of frequency response plot where the stable and saddletype bifurcation branches are crossed with each other, this arises due to the presence of directly excited first mode (principal parametric resonance) and indirectly excited second mode (3:1 internal resonance).

Conclusion
The current study aims to analyze the nonlinear transverse vibration of axially accelerating viscoelastic beam moving over elastically supported pulley mounting system under the joint influence of parametric as well as internal resonance. The axial speed is taken as harmonically varying component superimposed over a mean velocity while the longitudinal tension is the function of velocity, spatial coordinate, and pulley stiffness parameter. MMS method and subsequently continuation algorithm is implemented to study impact of various system parameters on the stability and bifurcation of the system.
The variation in pulley stiffness parameter brings significant qualitative and quantitative changes in nonlinear characteristics of the traveling system which is demonstrated in the form of stability and bifurcation features in the frequency and amplitude response plots. These results are unique, interesting and also presented here for the first time.
Apart from it the effect of variation of internal and external viscous damping on stability and bifurcation is also studied. Notable and interesting results are obtained which are illustrated in the previous two sections. Out of them some selected results are summarized below.
• Frequency and amplitude response plots display three types of equilibrium solution, i.e., trivial, nontrivial and nontrivial isolated solutions curves for higher values of the support stiffness parameter (k). However, for lower values of this parameter, typically when k ¼ 0:5 the third kind of solution curve disappears in the frequency response plot. • The maximum amplitude of response plots and the number of pitchfork bifurcation points remain unaffected with the decrease in k value, but the number of Hopf bifurcation points reduces along the nontrivial and nontrivial isolated curves. This results in the increase of Hope-type unstable portion in those solution curves. The same trend is observed for the trivial branch solution. Moreover, corresponding to the case of k ¼ 0:7; two additional saddle-type bifurcation points appear on the nontrivial branch (P 2 -P 3 ). • With increase in material damping a mixed impact is observed in the equilibrium solution of the frequency response analysis. Typically, when a increases to 0.0005 from 0, the number of saddle node bifurcation points increases to 4 from 3, but at another higher value of a = 0.001, it decreases to 2. The number of pitch-fork bifurcation points reduces from four to two in the last two cases. The range of Hopf-type instability in the equilibrium solution curves shrinks with the increase in a.
It may be noted that typically for a = 0.001, the nontrivial isolated steady state curve disappears in the frequency response plots. • The increase in external damping decreases the range of instability occurred due to Hopf bifurcation points. The number of pitchfork bifurcation points remains 4 for external damping (l) values 0.1 and 0.5, but when it changes to 0.7, the number decreases to 2. It may be noted that the impact of a on variation of nonlinear characteristics is more pronounced compared to that of l. • The system exhibits sensitiveness toward initial conditions while evaluating dynamic solutions. This is due to the multi-branch equilibrium solutions typical of the nonlinear system. Stable periodic, mixed mode, quasi-periodic and unstable chaotic behaviors are displayed by the traveling system at different initial points with a fixed set of control parameters. The above results based on parametric study may be useful for design and operation of traveling viscoelastic beams in various engineering and industrial applications.
Funding This research received no specific grant from any funding agency in the public, commercial, or not-for-profit sectors.

Declarations
Conflict of interest The authors declare that there is no conflict of interests regarding the publication of this paper.
Data availability All the data generated or analyzed during this study are included in this manuscript.