Internal resonance of an axially transporting beam with a two-frequency parametric excitation

This paper investigates the transverse 3:1 internal resonance of an axially transporting nonlinear viscoelastic Euler-Bernoulli beam with a two-frequency parametric excitation caused by a speed perturbation. The Kelvin-Voigt model is introduced to describe the viscoelastic characteristics of the axially transporting beam. The governing equation and the associated boundary conditions are obtained by Newton’s second law. The method of multiple scales is utilized to obtain the steady-state responses. The Routh-Hurwitz criterion is used to determine the stabilities and bifurcations of the steady-state responses. The effects of the material viscoelastic coefficient on the dynamics of the transporting beam are studied in detail by a series of numerical demonstrations. Interesting phenomena of the steady-state responses are revealed in the 3:1 internal resonance and two-frequency parametric excitation. The approximate analytical method is validated via a differential quadrature method.


Introduction
Axially transporting systems have wide applications in areas such as military, electronic, mechanical, and aerospace technology. The axial speed of a device may greatly affect the vibration characteristics of the transporting system, leading to variations of natural frequencies and complex modes. Therefore, it is of great significance to understand the vibration behavior * Citation: ZHANG, D. B., TANG, Y. Q., LIANG, R. Q., SONG, Y. M., and CHEN, L. Q. Internal resonance of an axially transporting beam with a two-frequency parametric excitation. Applied Mathematics and Mechanics (English Edition), 43 (12), 1805-1820 (2022) https://doi.org/10.1007/s10483-022-2930-9 of the transporting systems when they suffer from different excitation conditions, so that they may be safely employed in various engineering and industrial applications.
In the past decades, many researchers [1][2][3][4][5] have devoted themselves to investigating the transverse vibration of the axially transporting beam. Ding and Chen [6] developed the multiscale analysis to find the stability boundary of an axially accelerating beam, and introduced the finite difference schemes to verify the analytical results. Yang et al. [7] explored the effect of elastic foundation on the vibration and stability of an axially moving beam. Under a low and time-varying velocity, the transversal vibrations of a conveyor belt were discussed by Suweken and van Horssen [8][9][10] . Also, a similar study on an axially moving string with a time-varying velocity was done by Ponomareva and van Horssen [11] . By considering the effect of intermediate spring-support, the nonlinear stability and bifurcations of an axially accelerating beam were investigated by Ghayesh and Amabili [12] . Based on the Galerkin scheme and direct timeintegration, the three-dimensional nonlinear global dynamics of an axially moving viscoelastic beam was carried out by Farokhi et al. [13] . The kinematic aspects in modeling large-amplitude vibration of an axially moving beam were formulated in the mixed Eulerian-Lagrangian framework by Wang et al. [14] . Considering the effects of structural damping and axially moving, Ding et al. [15] studied the effects of the parameters on unstable regions and response curves. By considering the effects of finite deformation, Wang et al. [16] analyzed the vibration of an axially moving hyperelastic beam under the simply-supported condition.
As a result, in some practical cases, axially accelerating systems suffer from multi-frequency parametric excitations. In all of the above investigations, it was assumed that the axially transporting systems suffer from single-frequency excitations. The available papers [17][18][19][20] of an axially transporting system with a multi-frequency excitation are limited. The responses of twodegree-of-freedom and single-degree-of-freedom systems subject to a two-frequency parametric excitation were, respectively, determined by Nayfeh [21][22] . By considering the effect of speed fluctuations and multi-frequency tension, the parametric instability of axially moving media was investigated by Parker and Lin [23] . Three cases of an axially moving viscoelastic beam were studied with a multi-frequency excitation in Ref. [24]. By using the perturbation techniques, a nonlinear equation was obtained for a drillstring system in a deviated well with the timevarying axial load and axial velocity [25] . By considering the variable speed and axial force simultaneously, the dynamical stability of an axially moving beam was discussed byÖzhan [26] . A similar investigation was presented by Lv et al. [27] with the method of multiple scales to study the dynamic characteristics of an axially moving viscoelastic sandwich beam by considering the time-dependent axial tension and axially varying velocity.
All the above-mentioned investigations on the vibration of axially accelerating systems assumed the tension to be uniform along the longitudinal direction. In fact, the equal tension assumption is not accurate when the beam moves in acceleration. Chen and Tang [28][29] first proposed that the longitudinally varying tension due to the axial acceleration should be taken into account for axially accelerating systems. The nonlinear forced vibration of an axially accelerating viscoelastic tensioned beam with a two-frequency excitation was carried out by Ding and Zu [30] . The influence of the longitudinally varying tension on the nonlinear parametric vibration of an axially moving beam with a two-frequency parametric excitation has not been clear. To remedy the lack of investigation in this field, this paper aims to present the nonlinear vibration of an axially accelerating beam with a two-frequency parametric excitation and longitudinally varying tension.
Internal resonance is a typical nonlinear phenomenon. There are some important publications [31][32][33] on nonlinear transverse vibrations of an axially transporting system with internal resonance. The available papers of an axially transporting system with a multi-frequency excitation are limited. By considering the effect of the internal resonance and two-frequency parametric excitation, Sahoo et al. [34][35] studied the nonlinear transverse vibration of an axially accelerating viscoelastic beam with analytical and numerical methods. Nonlinear dynamics of an axially accelerating viscoelastic sandwich beam with a two-frequency parametric excitation and a 3:1 internal resonance were investigated by Zhu et al. [36] . Zhang et al. [37] studied the dynamic stability of an axially transporting beam with the internal resonance and two-frequency parametric excitation. So far, the nonlinear dynamic behavior of axial variable tension beam models with the internal resonance and two-frequency parametric excitation has not been understood.
In view of the lack of study on the axially transporting systems with the internal resonance and multi-frequency excitation, the present paper focuses on the steady-state responses of an axially accelerating viscoelastic beam subject to a 3:1 internal resonance and a two-frequency parametric excitation caused by speed perturbation. The study reveals that the nonlinear response of the transporting beam becomes more complex associated with the coupled effect of the 3:1 internal resonance and two-frequency parametric excitation.

Formulation
Consider a uniform horizontal axially transporting viscoelastic beam subject to a twofrequency parametric excitation in the presence of internal resonance. The model [37] of the beam with a time-dependent speed Γ(t), stretched along the x-axis to a length L by a longitudinally varying tension N x1 , is depicted in Fig. 1. The transverse displacement of any point on the beam at the neutral axis coordinate x is represented by the field variable v(x, t), where t is the time. The beam is modeled by the Euler-Bernoulli beam under simply-supported boundary conditions. For small-amplitude motions, it can be assumed that the transverse displacement and the longitudinal displacements are uncoupled. The governing equation of transverse vibration for the transporting beam is obtained from Newton's second law as where ρ is the density, A is the cross-sectional area, M (x, t) is the bending moment, k s is the stiffness of the foundation per unit length, and c d is the viscous damping. Considering the tension variation along the x-axis due to the axial acceleration, the longitudinally varying tension N x1 is obtained as where N x0 is the initial axial tension, and η is the axial support stiffness parameter varying between 0 (infinite rigidity) and 1 (no rigidity) [38] . η = 1/(1 + k r /(2EA/L)), in which k r is the relative stiffness of the support. The superscript dot denotes differentiation with respect to time. In the present investigation, the viscoelastic material of the beam obeys the Kelvin-Voigt model. Therefore, the moment-curvature relationship is expressed as where I is the area moment of inertia, E is Young's modulus, and α is the viscoelastic coefficient. For the axially transporting viscoelastic beam, the simply-supported boundary conditions are obtained as For the sake of generality, the non-dimensional variables and parameters are introduced as follows: where ε is a dimensionless parameter. It indicates that the transverse motion v, the viscoelasticity α, and the viscous damping c d are all small quantities. Then, substituting Eqs.
The non-dimensional inhomogeneous boundary conditions of the beam are obtained in the new parameters defined by Eq. (5) as The non-dimensional time dependent transport speed is expressed as where γ 0 denotes the non-dimensional mean speed, εγ n (n = 1, 2) denote the non-dimensional amplitudes, and Ω n (n = 1, 2) denote the non-dimensional frequencies of harmonically varying components. By substituting Eq. (8) into Eqs. (6) and (7), the nonlinear equation and the inhomogeneous boundary conditions of motion can be rewritten as The solutions to Eq. (9) can be expanded in the following form: where T 0 = t is the fast time scale, which presents vibrations with frequencies of the linear generating system. T 1 = εt is the slow time scale, which are used to obtain the responses and influence of the parametric excitation, the viscoelasticity of the material, and the nonlinearities. O(ε 2 ) denotes the higher-order terms of ε 2 , which are ignored in the following analysis. The time derivatives are By substituting Eqs. (11) and (12) into Eqs. (9) and (10), and equating the coefficients of ε with equal powers, we can obtain The solution to Eq. (13) may be assumed as where ϕ n (x) denotes the nth complex mode shape, and A n (T 1 ) presents an undetermined function. The nth complex frequency λ n includes two parts. One is the decay rate δ n , and the other is the frequency ω n . c.c. is the complex conjugate terms. Particularly, we can numerically solve the nth decay rate δ n and the natural frequency ω n from the following transcendental equation [37] : Meanwhile, the nth modal function can be obtained as With ξ = 1, κ = 0.5, and ζ = 1, Fig. 2 presents the variations of the first two natural frequencies of the beam with varying transporting velocities. It can be found from Fig. 2 that ω 2 is close to three times as much as ω 1 when γ 0 = 0.688 7 (ω 2 ≈ 3ω 1 , ω 1 = 3.225 355 49, and ω 2 = 9.676 049 25), resulting in a 3:1 internal resonance. Furthermore, we assume that there is no other commensurable frequency relationship with higher modes in this paper.

internal resonance
There are many possible choices for Ω 1 and Ω 2 resulting in resonances. In this paper, Ω 1 and Ω 2 are chosen to the first and second principal parametric resonances, respectively. For other choices of principal parametric resonance, the same technique can be used to determine stability conditions. Three frequency detuning parameters, i.e., σ 1 , σ 2 , and σ 3 , are used to quantify the proximity of ω 2 to 3ω 1 , Ω 1 to 2ω 1 , and Ω 2 to 2ω 2 . Then, the frequency relations can be expressed as By neglecting all terms not involved in resonances of Eq. (17), the solution to Eq. (13) can be truncated as Substituting Eqs. (20) and (21) into Eq. (15) yields where the overdot stands for the derivation of T 1 , the overbar indicates the complex conjugate, T ns denotes non-secular terms, and Internal resonance of an axially transporting beam with a two-frequency parametric excitation 1811 Generally speaking, the transverse vibration of gyroscopic continuum under homogeneous boundary conditions could be studied by applying the method of direct multiple scales. A general outline of the solvability conditions has been established by Chen and Zu [39] . Nevertheless, the boundary conditions become inhomogeneous, on account of the Kelvin viscoelastic constitutive relation. In fact, the problem for inhomogeneous boundaries can be solved theoretically by the approximate analytical method. Tang et al. [40] studied the nonlinear steady-state oscillating response of an axially transporting beam with internal resonance and inhomogeneous boundary conditions. In the present investigation, the same technique of the modified solvability conditions is presented to solve this problem. Then, the solution to Eq. (15) can be truncated to the first two modes, Substituting Eq. (24) into Eq. (22), and equalizing the coefficients of exp(iω 1 T 0 ) and exp(iω 2 T 0 ), we can obtain Using the distributive law, we have Based on the distribution law and the corresponding inhomogeneous boundary conditions (16), the left-hand sides of Eqs. (26a) and (26b) can be represented as Then, by applying the distributive law to Eqs. (27a) and (27b), we can obtain the two complex variable modulation equations as follows: According to the numerical results, ζ 1 and ξ 1 are positive real values. ζ 2 , ζ 3 , ζ 4 , ξ 2 , ξ 3 , and ξ 4 are complex numbers. ζ 5 , ζ 6 , ξ 5 , and ξ 6 are pure imaginary numbers with negative imaginary parts.
To investigate the equilibrium solutions and study the axially transporting system for stability and dynamic responses, the complex variable modulation equations (28a) and (28b) are introduced. The Cartesian transformation method is applied to denote the solutions to Eqs. (28a) and (28b) as in which p 1 , q 1 , p 2 , and q 2 are real functions of T 1 . S 1 = 0.5σ 2 , S 2 = 1.5σ 2 − σ 1 , and σ 3 = 3σ 2 − 2σ 1 . By substituting Eq. (29) into Eqs. (28a) and (28b), the normalized reduced equations are obtained aṡ Based on the normalized reduced equations (30a)-(30b), the nonlinear steady-state oscillating solutions and the dynamic behaviors of the axially transporting system can be numerically solved. In Eqs. (30a)-(30d),ṗ 1 =q 1 =ṗ 2 =q 2 = 0, when the system is stable. Namely, the terms of the left-hand side in Eqs. (30a)-(30d) are zero. Then, we can solve the nonlinear algebraic equations. In addition, the stabilities of the trivial solutions are obtained by the Lyapunov theory.

Analytical results and discussion
In the current research, the parameters are selected as α = 0.000 01, c d = 0.001, k N = 23.094 01, and γ 0 = 0.86. The related parameter values are obtained as σ 1 = 0.660 901 84, ω 11 = 2.917 115 43, and ω 12 = 9.412 248 13. The amplitudes of the harmonically varying components of the first two modes are γ 1 = 0.04 and γ 2 = 0.03. In the following investigation, only the positive sides of the frequency response curves are presented due to the fact that the frequency response curves are symmetric about the σ 2 -axis.
Typical frequency response curves and their enlarged view with the axial tension fluctuation frequency σ 2 for the first and second modes are presented in Fig. 3. The transporting beam system parameters are selected as α = 0.000 01 and c d = 0.001. In Fig. 3, the stable solutions are represented by solid lines, and the unstable solutions are displayed by dashed lines. At the same time, different colors are introduced to distinguish between different lines. As we can see in the figure, the transporting beam system displays a hardening-spring nonlinear feature. Figure 3 shows that only trivial solutions are present and are stable for the detuning parameter σ 2 < P 1 (σ 2 = −0.133 41). As σ 2 increases in the lower area, the trivial solution becomes unstable through a supercritical pitchfork bifurcation point P 1 , where a pair of complex conjugate eigenvalues passes through the imaginary axis of the complex plane from the left side to the right side. After that, trivial solutions lose their stability, and a stable nontrivial solution branch appears. With the increase in σ 2 , the stable nontrivial solution branch of the first mode from P 1 first increases and then decreases. As σ 2 increases, the unstable trivial solution regains stability through the pitchfork bifurcation point P 2 (σ 2 = 0.132 97). At the same time, another nontrivial solution branch appears at P 2 . With the increase in σ 2 from P 2 , the nontrivial solution branches of the first and second modes are always unstable. As σ 2 increases further, the trivial solutions lose stability again through another supercritical pitchfork bifurcation point at P 3 (σ 2 = 0.313 52). With the increase in σ 2 from P 3 , the stable nontrivial solution becomes unstable at a Hopf bifurcation point P 5 (σ 2 = 42 427). Then, the unstable trivial solutions become stable through a subcritical pitchfork bifurcation point at P 4 (σ 2 = 0.568 14). With the increase in σ 2 from P 4 , the nontrivial solution branch of the first mode is always unstable and increases monotonously. The nontrivial solution branch of the second mode from P 4 first increases and then decreases. When σ 2 1.044 31, an independent nontrivial solution branch (green line) appears. The nontrivial solution of the upper branch of green line becomes stable through a reverse Hopf bifurcation point P 6 (σ 2 = 1.126 46). After that, the nontrivial solution loses stability at a Hopf bifurcation point P 7 (σ 2 = 1.  Fig. 4(b), it can be easily found that the instability intervals are the same in the first and second modes. Inspecting Fig. 4, it can be found that the stable nontrivial solution branch at the supercritical pitchfork bifurcation point P 1 (P 3 ) and the unstable nontrivial solution branch at the subcritical pitchfork bifurcation point P 2 (P 4 ) form a solution curve, and its stability changes at a reverse Hopf bifurcation point P 5 (Hopf bifurcation point P 6 ). Different from Fig. 3, with the increase in the detuning parameter σ 2 , there are two independent nontrivial solution branches (blue and green curves), and their stability changes at P 7 and P 8 , respectively. In the first mode, the upper branch of the blue nontrivial solution curve is unstable, while the lower branch is stable. However, the upper branch of the green nontrivial solution curve is stable, and the lower branch is unstable. In the second mode, the variation trend is opposite. It can be seen from Fig. 5 that there are still four pitchfork bifurcation points containing P 1 , P 2 , P 3 , and P 4 on the trivial solution branch. The stable nontrivial solution branch at the supercritical pitchfork bifurcation point P 1 (P 3 ) and the unstable nontrivial solution branch at the subcritical pitchfork bifurcation point P 2 (P 4 ) form a solution curve, and its stability changes at the reverse Hopf bifurcation point P 5 (Hopf bifurcation point P 6 ). With the increase in σ 2 , two independent nontrivial solution branches (blue and green curves) appear, and their stability changes at the saddle node bifurcation point P 7 and the reverse Hopf bifurcation point P 8 , respectively. Different from Fig. 4, the blue nontrivial solution curve is above the green nontrivial solution curve with the increase in σ 2 .
Comparing these response curves shown in Fig. 5, there are significant changes for the response curves in Fig. 6. It can be easily found from Fig. 6 that the stable nontrivial solution branch at the supercritical pitchfork bifurcation point P 1 loses stability through the Hopf bi- furcation point P 5 . Then, it regains stability through a reverse Hopf bifurcation point P 7 . The nontrivial solution branch at the subcritical pitchfork bifurcation point P 2 is always unstable. The stable nontrivial solution branch at P 3 and the unstable nontrivial solution branch at P 4 form a solution curve, and its stability changes at a saddle node bifurcation point P 6 . With the increase in σ 2 , an independent nontrivial solution branch (red curve) appears, and its stability changes at the saddle node bifurcation point P 8 . In the first mode, the independent nontrivial solution branch decreases monotonously with σ 2 increasing. However, in the second mode, the variation trend is opposite. Comparing with the response curves presented in Fig. 6, the larger instability frequency interval disappears with the increase in the viscoelastic coefficient α in Fig. 7. In addition, there is an isolated loop in the local region of the non-zero solution. The stability of the nontrivial solution curve changes at the saddle node bifurcation point P 5 and Hopf bifurcation point P 6 , respectively. With the increase in the viscoelastic coefficient α, the isolated loop disappears in Fig. 8. The amplitude of the nontrivial branch of the first mode increases monotonically in Fig. 8(a). However, the amplitude of the second mode has a peak value in Fig. 8(b). The unstable amplitudes of the second mode increase first and then decrease. The variation of the stable amplitude is more complex in the second mode. Comparing Figs the material has more influence on the steady-state vibration responses. The vibration energy of the transporting system will continue to be dissipated due to the introduction of Kelvin's viscoelastic constitutive relation. The stability of the transporting system during long-term vibration is improved based on the irreversibility of viscous deformation.

Numerical validation
In this section, the same form of sampling points with Tang and Zhang [40] is selected. The differential quadrature analogue of the governing equation (6) can be obtained as ik − (x i − 1)(γ 1 Ω 1 cos(Ω 1 t) + γ 2 Ω 2 cos(Ω 2 t)) A With v 0 = 0.004 sin(πx i ), σ 2 = 0.418, N = 15, c d = 0.001, α = 0.000 01, γ 0 = 0.86, γ 1 = 0.04, and γ 2 = 0.03, Fig. 9 presents the time histories and phase portraits of the beam center (x i = 0.5). Figure 10 shows partial enlarged views of Fig. 9. The amplitude spectrum obtained by the the fast Fourier transform (FFT) with t 1 = 4 000, t 2 = 5 000, and T = 0.001 is shown in Fig. 11. Inspecting Fig. 11, there are two special frequency points ω 3 = 3.126 and ω 4 = 9.379. The amplitudes of the first and second modes are much larger than those of the other modes. Namely, the first-and second-order modes play a major role at this moment. To contrast the quantitative difference between the numerical results with the approximate analytical outcomes for the response amplitudes of the midpoint of the axially transporting beam, the amplitude and phase of every mode must be considered. In the current investigation, with N = 15, α = 0.000 01, c d = 0.001, γ 0 = 0.86, γ 1 = 0.04, and γ 2 = 0.03, the comparison for the analytical and numerical outcomes of the response amplitudes is presented in Fig. 12. It can be easily found that the response amplitudes have complex dynamic behaviors. In addition, the comparison illustrates that the outcomes of the two methods are in good agreement.

Conclusions
This paper treats nonlinear transverse vibrations of an axially transporting viscoelastic beam in the presence of the 3:1 internal resonance and two-frequency parametric excitation. The governing equation and the associated boundary conditions are obtained from Newton's second law. The method of multiple scales is developed to study the steady-state responses along with their stabilities and bifurcations. The numerical examples are presented to illustrate the effect of the viscoelastic coefficients of the material on the dynamic behaviors. The validity of the analytical results is demonstrated by the numerical results via a differential quadrature method. The main conclusions are listed as follows.
(i) There are two types of steady-state responses, i.e., trivial zero solutions and nontrivial solutions. The transporting system displays a hardening-spring nonlinear feature. Duo to the 3:1 internal resonance, energy exchange is discovered by both analytical and numerical analysis. These phenomena are due to modal interactions caused by the internal resonance. Even at the primary resonance, nonlinear modal interactions result in an important variation of the steadystate responses.
(ii) The complex frequency response curves are detected with the co-existence of the 3:1 internal resonance and two-frequency parametric excitation. The instability frequency interval decreases with the increase in the material viscoelastic coefficient. In addition, the supercritical and subcritical pitchfork bifurcation points of the trivial solution disappear gradually. The vibration energy of the transporting system will continue to be dissipated due to the viscoelasticity constituted by Kelvin's relation. The stabilities of the transporting system during long-term vibration are improved based on the irreversibility of viscous deformation.
(iii) The approximate analytical stable steady-state responses for principal parametric resonances of the midpoint of the transporting beam represent good agreement in qualitatively and quantitatively with those obtained by the numerical integrations.