Nonlinear Vibrations of Axially Moving Beams

1. Introduction

Axially moving beams can represent many engineering devices, such as band saws, power transmission belts, aerial cable tramways, crane hoist cables, flexible robotic manipulators, and spacecraft deploying appendages. Despite usefulness and advantages of these devices, vibrations associated with the devices have limited their applications. Therefore, understanding transverse vibrations of axially moving beams is important for the design of the devices. The investigations on vibrations of axially moving beams have theoretical importance as well, because an axially moving beam is a typical representative of distributed gyroscopic systems. The term “gyroscopic” arose in recognition of an early problem in gyrodynamics. Actually, the Coriolis acceleration component experienced by axially moving materials imparts a skew-symmetric or gyroscopic term to their governing equations. Due to particular characteristics of the gyroscopic term, the approaches developed in analyzing vibrations of an axially moving string can be applied to other more complicated distributed gyroscopic systems. Because of the practical and theoretical significance, the investigation on nonlinear vibrations of axially moving beams is a challenging subject which has been studied for many years and is still of interest today.

The relevant researches on transverse vibrations of axially moving strings can be dated back to (Aiken, 1878). There are several excellent and comprehensive survey papers, notably (Mote, 1972), (Ulsoy and Mote, 1978), (Mote et al., 1982), (D’Angelo et al., 1985), (Wickert and Mote, 1988), (Wang and Liu, 1991), (Abrate, 1992), (Zhu, 2000), reviewing the state-of-the-art in different time phases of investigations related to vibrations of axially moving beams. The present chapter emphasizes on the recent achievements, although some early results are mentioned for the sake of completeness. Besides, the chapter focuses the nonlinear problem only. If the vibration amplitude is large, the nonlinearity should be taken into account. Hence the chapter, unlike (Chen, 2005a) for axially moving strings, is not a comprehensive survey with a complete and detailed representation of current researches. Instead, the chapter is a counterpart of (Chen et al., 2008) for axially moving beams. The author tries to put the some available results into a general framework, as well as to highlight the work of the author and his collaborators. It is hoped that the chapter serves as a collection of ideas, approaches, and main results in investigations on nonlinear vibration of axially moving beams.

The chapter is organized as follows: Section 2 focuses on the mathematical models of nonlinear transverse vibration. The special attentions are paid to the comparison of two different nonlinear models and the introduction of the material time derivative into the viscoelastic constitutive relations. Section 3 covers the developments and the applications of approximately analytical methods, including the asymptotic method, the Lindstedt-Poincaré method, the method of normal forms, the method of nonlinear, complex modes, the method of multiple scales, and the incremental harmonic balance method. Section 4 is devoted to the numerical approaches, including the Galerkin discretization, the finite difference, and the differential quadrature. Section 5 reveals the nonlinear behaviors such as bifurcation and chaos based on the numerical solutions. Section 6 discusses energetics, conserved quantity and the applications. The final section recommends future research directions.

2. Governing equations

2.1. Coupled vibration

The governing equation is the base of all analytical or numerical investigations. Generally, an axially moving beam undergoes both the longitudinal vibration and the transverse vibration, and they are coupled. (Thurman & Mote, 1969) obtained the governing equation of coupled longitudinal and transverse vibrations of an axially moving beam. (Koivurova & Salonen 1999) revisited the same modeling problem and clarified its kinematic aspects. Their nonlinear formulation for the moving beam problem has two limitations: the material of the beam is linear elastic constituted by Hooke’s law, and the axial speed of the beam is a constant. As (Wickert & Mote 1988) pointed out, modeling of dissipative mechanisms is an important vibrations analysis topic of axially moving materials. An effective approach is to model the beam as a viscoelastic material. Therefore, it is necessary to deal with constitutive laws other than Hooke’s law. Axial transport acceleration frequently appears in engineering systems. For example, if an axially moving beam models a belt on a pair of rotating pulleys, the rotation vibration of the pulleys will result in a small fluctuation in the axial speed of the belt. The nonlinear model in (Thurman & Mote, 1969) for coupled vibration can be generalized to an axially accelerating viscoelastic beam as follows.

Consider a uniform axially moving beam of density ρ, cross-sectional area A, moment of inertial I, and initial tension P₀, as shown in Figure 1. The beam travels at the uniform transport speed γ between two boundaries separated by distance l. Assume that the deformation of the beam is confined to the vertical plane. A mixed Eulerian-Lagrangian description is adopted. The distance from the left boundary is measured by fixed axial coordinate x. The beam is subjected to an external excitation f_u(x,t) and f_v(x,t) in longitudinal and transverse directions respectively, where t is the time. The in-plane motion of the beam is specified by the longitudinal displacement u(x,t) related to coordinate translating at speed γ and the transverse displacement v(x,t) related to the spatial frame.

Figure 1.

The physical model of an axially accelerating beam.

Study the motion of the beam in a reference frame moving in the axial direction and at speed γ. The reference system is not an inertial frame if γ is not a constant. The beam is a one-dimensional continuum undergoing an in-plane motion in the moving reference frame, the Eulerian equation of motion of a continuum gives

$$\begin{array}{l}\rho \frac{d^{2}u}{d{t}^2}=\frac{\partial }{\partial x}\left(\frac{1}{A}\frac{\left(P_0+A\sigma \right)\left(1+u{,}_x\right)}{\sqrt{{\left(1+u{,}_x\right)}^2+v{,}_x^2}}\right)-\rho \dot{\gamma }+\frac{f_u\left(x,t\right)}{A},\\ \rho \frac{d^{2}v}{d{t}^2}=\frac{\partial }{\partial x}\left(\frac{1}{A}\frac{\left(P_0+A\sigma \right)v{,}_x}{\sqrt{{\left(1+u{,}_x\right)}^2+v{,}_x^2}}\right)-\frac{M{,}_{xx}\left(x,t\right)}{A}+\frac{f_v\left(x,t\right)}{A}\end{array}$$E1

,

where a comma preceding x or t denotes partial differentiation with respect to x or t, σ(x,t) is the axial disturbed stress, and M(x,t) is the bending moment. The viscoelastic material of the beam obeys the Kelvin model, with the constitution relation

$$\sigma \left(x,t\right)=\left(E+\eta \frac{d}{dt}\right){\epsilon }_N\left(x,t\right)$$E2

,

where, E is the Young's modulus, η is the dynamic viscosity, and the disturbed strain ε_N(x, t) of the beam is given by the nonlinear geometric relation

$${\epsilon }_N=\sqrt{{\left(1+u{,}_x\right)}^2+v{,}_x^2}-1$$E3

For a slender beam (for example, with I/(Al²)<0.001), the linear moment-curvature relationship of Euler-Bernoulli beams is sufficiently accurate,

$$M\left(x,t\right)=\left(E+\eta \frac{d}{dt}\right)Iv{,}_{xx}$$E4

In the moving reference frame, the beam itself is without any axial transportation, while the boundaries are moving at speed γ. The axially moving beam is constrained by rotating sleeves with rotational springs (Chen & Yang, 2006a). The stiffness constant of two springs is the same, denoted as K. Nullifying the transverse displacements and balancing the bending moment at both ends lead to the boundary conditions

$$u\left(s,t\right)=0,u\left(l+s,t\right)=0$$E5

,

$$v\left(s,t\right)=0,EIv{,}_{xx}\left(s,t\right)-Kv{,}_x\left(s,t\right)=0;v\left(l+s,t\right)=0,EIv{,}_{xx}\left(l+s,t\right)+Kv{,}_x\left(l+s,t\right)=0$$E6

.

where $\dot{s}=\gamma$. To avoid the moving boundary conditions (5), which are difficult to tackle, the transformation of coordinates is introduced as follows

$$x\leftrightarrow x+s,t\leftrightarrow t$$E7

.

Then, expressed in the new coordinates, the boundary conditions have a simpler form

$$u\left(0,t\right)=0,u\left(l,t\right)=0$$E8

;

$$v\left(0,t\right)=0,EIv{,}_{xx}\left(0,t\right)-Kv{,}_x\left(0,t\right)=0;v\left(l,t\right)=0,EIv{,}_{xx}\left(l,t\right)+Kv{,}_x\left(l,t\right)=0$$E9

.

Under the new coordinates, the partial derivatives with respect to x and t remain invariant, and the total time derivative changes as follows

$$\frac{d}{dt}\leftrightarrow \gamma \frac{\partial }{\partial x}+\frac{\partial }{\partial t}$$E10

Substitution of equations (2), (4), and (10) into equation (1) yields

$$\begin{array}{l}\rho A\left(u{,}_{tt}+2\gamma u{,}_{xt}+\dot{\gamma }\left(1+u{,}_x\right)+{\gamma }^{2}u{,}_{xx}\right)=\\ \frac{\partial }{\partial x}\left(\frac{\left[P_0+A\left(E{\epsilon }_N+\eta {\epsilon }_N{,}_t+\eta \gamma {\epsilon }_N{,}_x\right)\right]\left(1+u{,}_x\right)}{\sqrt{{\left(1+u{,}_x\right)}^2+v{,}_x^2}}\right)+f_u\left(x,t\right),\\ \rho A\left(v{,}_{tt}+2\gamma v{,}_{xt}+\dot{\gamma }v{,}_x+{\gamma }^{2}v{,}_{xx}\right)=\\ \frac{\partial }{\partial x}\left(\frac{\left[P_0+A\left(E{\epsilon }_N+\eta {\epsilon }_N{,}_t+\eta \gamma {\epsilon }_N{,}_x\right)\right]v{,}_x}{\sqrt{{\left(1+u{,}_x\right)}^2+v{,}_x^2}}\right)-\left[EIv{,}_{xxxx}+\eta I\left(v{,}_{xxxxt}+\gamma v{,}_{xxxxx}\right)\right]+f_v\left(x,t\right)\end{array}$$E11

,

If other viscoelastic constitutive relations are used to describe the beam materials, they can be incorporated into the governing equation in the similar way. However, a controversial issue arises concerning the application of differential-type constitutive laws including the Kelvin relation in axially moving materials. Some investigators used the partial time derivative in the Kelvin model for axially moving strings (Zhang & Zu, 1998) (Zhang & Song, 2007), (Chen et al., 2007) and (Ghayesh, 2008), or beams (Chen & Yang, 2005, 2006a,b), (Ghayesh & Balar, 2008), (Ghayesh & Khadem, 2008), (Yang et al., 2009), and (Özhan & Pakdemirli, 2009). However, (Mochensturm & Guo, 2005) convincingly argued that the Kelvin model generalized to axially moving materials should contain the material time derivative to account for the added “steady state” dissipation of an axially moving viscoelastic string. Actually the material time derivative was also used in other works on axially moving viscoelastic beams (Marynowski, 2002, 2004, 2006), (Marynowski & Kapitaniak, 2002, 2007), (Yang & Chen, 2005), (Ding & Chen, 2008), (Chen & Ding, 2008, 2009), (Chen & Wang, 2009) and (Chen, et al., 2008, 2009, 2010). Here a coordinate transform will be proposed to develop the governing equations, which can introduce naturally the material time derivative in the viscoelastic constitutive relations.

In small but finite stretching problems in literatures of nonlinear oscillations, only the lowest order nonlinear terms need to be retained so that the governing equation of small-amplitude motion will be obtained. Such simplified coupled governing equations were used in analytical investigations on axially moving elastic beams (Thurman & Mote, 1969), (Riedel & Tan, 2002), and (Sze et al., 2005).

It should be remarked that there are different types of governing equations for axially moving beams (Tabarrok et al., 1974), (Wang & Mote, 1986, 1987), (Wang, 1991), (Hwang & Perkins, 1992a,b, 1994), (Vu-Quoc & Li 1995), (Behdinan, et al, 1997), (Hochlenert et al., 2007), (Pratiher & Dwivedy 2008), (Spelsbrg-Korspeter et al., 2008), and (Humer & Irschik, 2009). Actually, there are various beam theories such as Euler-Bernoulli theory, shear-deformable theories, and three-dimensional theories, and geometric nonlinearities may take different forms. Correspondingly, there are various governing equations of axially moving beams. Even if an axially stationary slender structure is prescribed by more sophisticated governing equations, the coordinate transform (7) is a still convenient approach to derive the governing equations of the slender structure undergoing an axial motion.

2.2. Transverse vibration

Although the transverse vibration is generally coupled with the longitudinal vibration, many researchers considered only the transverse vibration in order to derive a tractable equation. Inserting u=0 into equation (3) and then omitting higher order nonlinear terms yield a simplified strain-displacement relation termed as the Lagrange strain

$${\epsilon }_L=v{,}_x^2/2$$E12

Inserting u=0 into equation (11) and then retaining lower order nonlinear terms only yield a nonlinear partial-differential equation

$$\begin{array}{l}\rho A\left(v{,}_{tt}+2\gamma v{,}_{xt}+\dot{\gamma }v{,}_x+{\gamma }^{2}v{,}_{xx}\right)-P_{0}v{,}_{xx}+\left[EIv{,}_{xxxx}+\eta I\left(v{,}_{xxxxt}+\gamma v{,}_{xxxxx}\right)\right]\\ =\frac{\partial }{\partial x}\left[\left(AE{\epsilon }_L+A\eta {\epsilon }_L{,}_t+A\eta \gamma {\epsilon }_L{,}_x\right)v{,}_x\right]+f_v\left(x,t\right)\end{array}$$E13

.

The quasi-static stretch assumption means that one can use the averaged value of the disturbed tension $${\displaystyle {\int }_0^1\left[AE{\epsilon }_L+A\eta \left({\epsilon }_L{,}_t+\gamma {\epsilon }_L{,}_x\right)\right]dx}/l$$ to replace the exact value AEε+Aη(ε,_t+cε,_x). Thus equation (18) leads to nonlinear integro-partial-differential equation

$$\begin{array}{l}\rho A\left(v{,}_{tt}+2\gamma v{,}_{xt}+\dot{\gamma }v{,}_x+{\gamma }^{2}v{,}_{xx}\right)-P_{0}v{,}_{xx}+\left[EIv{,}_{xxxx}+\eta I\left(v{,}_{xxxxt}+\gamma v{,}_{xxxxx}\right)\right]\\ =\frac{v{,}_{xx}}{l}{\displaystyle {\int }_0^l\left(AE{\epsilon }_L+A\eta {\epsilon }_L{,}_t+A\eta \gamma {\epsilon }_L{,}_x\right)}dx+f_v\left(x,t\right)\end{array}$$E14

.

Both equation (13) and equation (14) are governing equations of transverse nonlinear vibration.

Both the nonlinear partial-differential equation and the nonlinear integro-partial-differential equation have been applied to some special cases such as free vibration without external excitation (F_v=0), elastic beams without viscoelasticity (η=0), uniformly moving beams without axially acceleration ($\dot{\gamma }$=0). The applications of the nonlinear partial-differential equation include (Chen & Zu, 2004) for uniformly moving elastic beams without external excitation, (Marynowski, 2002, 2004) and (Marynowski & Kapitaniak, 2007) for axially moving viscoelastic beams without external excitation, (Yang & Chen, 2005) and (Chen & Yang, 2006) for axially accelerating viscoelastic beams, and (Chen et al., 2007) for uniformly moving elastic beams without external excitation. The applications of the nonlinear integro-partial-differential equation include (Wickert, 1991), (Pellicano & Zirilli, F., 1997), (Pellicano & Vestroni, 2000), (Chakraborty & Mallik, 2000a), (Pellicano, 2001), (Kong & Parker, 2004) and (Chen & Zhao, 2005) for uniformly moving elastic beams without external excitation, (Ghayesh, 2008) for uniformly moving viscoelastic beams without external excitation, (Pellicano & Vestroni, 2000), (Özhan & Pakdemirli, 2009) for uniformly moving elastic beams, (Chakraborty & Mallik, 1999), (Öz et al, 2001) and (Ravindra & Zhu, 1998) for axially accelerating elastic beams without external excitation, (Chakraborty & Mallik, 1998) (Chakraborty et al., 1999), (Chakraborty & Mallik, 2000b) for axially moving elastic beams, (Parker & Lin, 2001) for axially accelerating elastic beams, and (Yang et al., 2009), (Chen et al., 2009) for axially accelerating viscoelastic beams, and (Özhan & Pakdemirli, 2009) for uniformly moving viscoelastic beams. Approximately analytical investigations on free vibration of axially moving elastic (Chen & Yang, 2007), forced vibration of axially moving viscoelastic beams (Yang & Chen, 2006), and parametric vibration of axially accelerating viscoelastic beams (Chen & Yang, 2005) and (Chen & Ding, 2008) demonstrated that the nonlinear partial-differential equation and the nonlinear integro-partial-differential equation yield the qualitatively same results but there are quantitative differences.

The nonlinear integro-partial-differential equation can also be obtained through uncoupling the governing equation for coupled longitudinal and transverse vibration under the quasi-static stretch assumption in small but finite stretching problems, and a special case of free vibration of axially moving elastic beam was treated in (Wickert, 1992). Under quasi-static stretch assumption, the dynamic tension to be a function of time alone. In traditional derivation in (Wickert, 1992), the nonlinear integro-partial-differential equation seems more exact than the nonlinear partial-differential equation because it is the transverse equation of motion in which the longitudinal displacement field is taken into account. However, the derivation here indicates that the nonlinear partial-differential equation can be reduced to the nonlinear integro-partial-differential equation based on the quasi-static stretch assumption. Numerical investigations on free vibration of axially moving elastic beams (Ding & Chen, 2008) and forced vibration of axially moving viscoelastic beams (Chen & Ding, 2009) indicated that the nonlinear integro-partial-differential equation is superior to the partial-differential equation, in the sense that approximates the coupled governing equation of planar motion better (some details in Subsection 4.2). However, since there has no decisive evidence to favor any models of transverse nonlinear vibration of axially moving beams, it is still an open problem.

3. Approximate analytical methods

3.1. Direct-perturbation approaches

As exact solutions are usually unavailable, approximate analytical methods are widely applied to investigate nonlinear vibration of axially moving beams. The approximate analytical methods can be applied to the nonlinear (integro-)partial-differential equations without discretization. Such a treatment is regarded as a direct-perturbation. The practice can be dated back to (Thurman & Mote, 1969) in which a modified Lindstedt method was used to calculate the fundamental frequency.

The method of multiple scales can be employed to analyze nonlinear vibration of axially moving beams. Actually, a general framework of the multi-scale analysis has been proposed for a linear gyroscopic continuous system under small nonlinear time-dependent disturbances (Chen & Zu, 2008). Consider a gyroscopic continuous system with a weak disturbance

$$Mv{,}_{tt}+Gv{,}_t+Kv=\epsilon N\left(x,t\right)$$E15

.

where v(x,t) is the generalized displacement of the system at spatial coordinate x and time t, linear, time-independent, spatial differential operators M, G and K represent mass, gyroscopic and stiffness operators respectively, ε stands for a small dimensionless parameter, and N(x,t) expresses a nonlinear function of x and t that may explicitly contain v and its spatial and temporal partial derivatives as well as its integral over a spatial region or a temporal interval. N(x,t) is periodic in time with the period 2π/ω. Define an inner product

$$⟨f,g⟩={\displaystyle \underset{E}{\int }f\left(x\right)\overline{g}\left(x\right)dx}$$E16

,

for complex functions f and g defined in the gyroscopic continuum E, where the overbar denotes the complex conjugate. M, K are symmetric and G is skew symmetric in the sense

$$⟨Mf,g⟩=⟨f,Mg⟩,⟨Kf,g⟩=⟨f,Kg⟩,⟨Gf,g⟩=-⟨f,Gg⟩$$E17

for all functions f and g satisfying appropriate boundary conditions. A uniform approximation is sought in the form

$$v\left(x,t\right)=v_0\left(x,T_0,T_1\right)+\epsilon v_1\left(x,T_0,T_1\right)+O\left({\epsilon }^2\right)$$E18

,

where T₀=t, T₁=εt, and O(ε²) denotes the term with the same order as ε² or higher. Substitution of equation (18) into equation (15) yields

$$M{v}_0{,}_{T_0{T}_0}+G{v}_0{,}_{T_0}+K{v}_0=0$$E19

,

$$M{v}_1{,}_{T_0{T}_0}+G{v}_1{,}_{T_0}+K{v}_1=N_1\left(x,T_0,T_1\right)$$E20

,

where N₁(x,T₀,T₁) stands for a nonlinear function of x, T₀ and T₁, which usually depends explicitly on v₀ and its derivatives and integrals. In addition, N₁(x,T₀,T₁) is periodic in T₀ with the period 2π/ω. Separation of variables leads to the solution of equation (19) as

$$v_0\left(x,T_0,T_1\right)={\displaystyle \sum _{j=1}^{\infty }{A}_j\left(T_1\right){\varphi }_j\left(x\right)e^{i{\omega }_j{T}_0}}+cc$$E21

,

where A_j denotes a complex function to be determined later, φ_j and ω_j represents respectively the complex modal function and the natural frequency given by

$$-{\omega }_j{}^{2}M{\varphi }_j+i{\omega }_{j}G{\phi }_j+K{\phi }_j=0$$E22

and the boundary conditions, and cc stands for the complex conjugate of all preceding terms on the right side of the equation. If ω approaches a linear combination of natural frequencies of equation (19), the summation parametric resonance may occur. A detuning parameter σ is introduced to quantify the deviation of ω from the combination, and ω is described by

$$\omega ={\displaystyle \sum _{j=1}^{\infty }{c}_j{\omega }_j}+\epsilon \sigma$$E23

,

where c_j are real constants that are not all zero and only a finite of them are not zero. To investigate the summation parametric resonance, substitution of equations (21) and (23) into equation (20) leads to

$$M{v}_1{,}_{T_0{T}_0}+G{v}_1{,}_{T_0}+K{v}_1={\displaystyle \sum _{j=1}^{\infty }{F}_j\left(x,T_1\right)e^{i{\omega }_j{T}_0}}+NST+cc$$E24

,

where F_j(x, T₁) (j=1,2,…) are complex functions dependent explicitly on A_j(T₁) and their temporal derivatives as well as φ_j(x) and their spatial derivatives and integrals. (Chen & Zu, 2008) proved that the solvability condition is the orthogonality of the coefficient of the resonant term in the first order equation and the corresponding modal function of the zero order equation, e.g.

$$⟨F_j\left(x,T_1\right),{\phi }_j⟩=0$$E25

.

It should be noticed that the solvability condition (25) holds providing the boundary conditions are appropriate. That is, M and K are symmetric and G is skew symmetric under the boundary conditions. In a specific problem, these requirements can be checked for a given the operators, boundary conditions and the modal functions. However, the examination depends only on the unperturbed linear part of the problem. For example, equation (25) holds for an axially moving beam under condition (9) (Chen & Zu, 2008). Usually, it is assumed that only the modes involved in the resonance (23) need to be considered in the linear solution (21), and the assumption is physically sound. Some case studies demonstrated mathematically that the mode uninvolved in the resonance has no effect on the steady-state response (Ding & Chen, 2008), (Chen & Wang, 2009), and (Chen et al., 2009). (Özhan & Pakdemirli, 2009) proposed multi-scale analysis on forced vibrations of general continuous systems with cubic nonlinearities in the primary resonance case.

The method of multiple scales has been applied in various transverse nonlinear vibration problems of axially moving beams. These problems include free (Öz et al, 2001) and (Chen & Yang, 2007), forced (Özhan & Pakdemirli, 2009), and parametric(Öz et al, 2001) and (Özhan & Pakdemirli, 2009) vibration of axially moving elastic beams, as well as forced (Yang & Chen, 2006) and parametric (Chakraborty & Mallik, 1999), (Chen & Yang, 2005) and (Chen & Ding, 2008) vibrations of axially moving viscoelastic beams. In addition to these works on the base of the Euler-Bernoulli beam theory, the method of multiple scales has also be applied to study free vibration of an axially moving beam with rotary inertia and temperature variation effects (Ghayesh & Khadem, 2008), parametric vibration of axially moving viscoelastic Rayleigh beams (Ghayesh & Balar, 2008), and forced (Tang et al. 2009) and parametric (Tang et al., 2010) vibrations of axially moving elastic Timoshenko beams, while the multi-scale analysis on axially moving viscoelastic Timoshenko beams has been only limited to linear parametric vibration (Chen et al., 2010).

Addition to the method of multiple scales, the method of asymptotic analysis is also an effective approach to treat parametric or nonlinear vibration. Based on the idea of Krylov, Bogoliubov, and Mitropolsky, (Wickert, 1992) developed a asymptotic method for general gyroscopic continuous systems with weak nonlinearities, and the method was specialized to free nonlinear vibration of an axially moving elastic beam with supercritical transport speed. (Maccari, 1999) proposed another asymptotic approach for analyzing transverse vibration of axially stationary beams, which are disturbed conservative continuous systems, and determined external force-response and frequency-response curves in the cases of primary resonance and subharmonic resonance for a weakly periodically forced beam with quadratic and cubic nonlinearities. The approach was extended to the gyroscopic continuous system with a weak nonlinear and time-dependent disturbance in order to analyze transverse vibration of an axially accelerating viscoelastic string constituted by the Kelvin model (Chen et al., 2008) and the standard linear solid model (Chen & Chen, 2009). The method of asymptotic analysis has been also presented for nonlinear parametric vibration of axially accelerating viscoelastic beams constituted by the Kelvin model (Chen et al., 2009) as well as linear parametric vibration of axially accelerating viscoelastic beams constituted by the Kelvin model (Chen & Wang, 2009) and the standard linear solid model (Wang & Chen, 2010).

Nonlinear normal modes whose shapes depend on the amplitude provide a possible direct treatment on nonlinear vibration of axially moving beams. (Chakraborty et al., 1999) used a temporal harmonic balance and a spatial perturbation technique to determine the nonlinear complex normal modes for free and forced vibrations of axially moving elastic beams. The approach was adopted to study the response of a parametrically excited axially moving beam both without and with an external harmonic force (Chakraborty & Mallik, 1998). The results were justified by the wave propagation analysis (Chakraborty & Mallik, 2000a,b).

3.2. Discretization-perturbation approaches

Discretization of governing equations is a commonly used approach to obtain approximate solutions of vibration problems of continuous systems. For the governing equations (11) of coupled vibration of axially moving beams, one assumes an approximate solution in the form

$$u\left(x,t\right)={\displaystyle \sum _{i=1}^m{p}_i\left(t\right){\varphi }_i\left(x\right)}$$E26

,

$$v\left(x,t\right)={\displaystyle \sum _{i=1}^n{q}_i\left(t\right){\phi }_i\left(x\right)}$$E27

,

where p_i(t) and q_i(t) are generalized coordinates, and x _i(x) are φ _i(x) base functions that are usually chosen to be the linear vibration mode shapes of axially stationary beams or moving beams (Wickert & Mote, 1990), (Chen & Yang, 2006), and (Tang et al. 2008). A weighted-residual procedure such as the Galerkin procedure can be applied to truncate equation (11) into m+n nonlinearly coupled second-order ordinary-differential equations. A general description of the Galerkin procedure is as follows. Denote the differences between the left and right sides of two equations in equation (11) as F_u(x,u,v,t) and F_v(x,u,v,t), which are nonlinear functions of x and t that may explicitly contain v and its spatial and temporal partial derivatives as well as its integral over a spatial region or a temporal interval. Then approximate solution (31) and (32) satisfies

$$\begin{array}{l}⟨F_u\left(x,{\displaystyle \sum _{i=1}^m{p}_i\left(t\right){\varphi }_i\left(x\right)},v{\displaystyle \sum _{i=1}^n{q}_i\left(t\right){\phi }_i\left(x\right)},t\right),{\zeta }_j\left(x\right)⟩\\ =0,⟨F_v\left(x,{\displaystyle \sum _{i=1}^m{p}_i\left(t\right){\varphi }_i\left(x\right)},v{\displaystyle \sum _{i=1}^n{q}_i\left(t\right){\phi }_i\left(x\right)},t\right),{\psi }_k\left(x\right)⟩=0,\left(j=1,\cdots ,m;k=1,\cdots ,n\right)\end{array}$$E28

where ζ_j(x) and ψ_k(x) are the weight functions.

After the discretization, various perturbation techniques such as the method of multiple scales can be employed to analyze the resulting nonlinear ordinary-differential equations approximately. Such a treatment is regarded as a discretization-perturbation

In practical problems, m and n in the discretization expressions (26) and (27) are rather small, and they are usually 1 or 2. (Riedel & Tan, 2002) applied the method of multiple scales to the discretized equations (m=n=2) to determine the forced response of an axially moving elastic beam with internal resonance. The method of multiple scales was also applied to the discretized problem of coupled vibration of axially moving beams (Feng & Hu, 2002, 2003). (Sze et al., 2005) presented a general description of discretization of the governing equation of an axially moving elastic beam, and used incremental harmonic balance method to a concrete case of m=n=2 for forced response with internal resonance. In both studies, the mode shapes of axially stationary beams were chosen as the base functions and the weight functions.

Discretization-perturbation approaches have also been used in analyzing transverse nonlinear vibration of axially moving beams. In this case, equation (27) will be substituted into equation (13) or (14) and then the Galerkin procedure can be used to discretize equation (13) or (14) into n nonlinearly coupled second-order ordinary-differential equations that can be solved approximately via various perturbation techniques. The Lindstedt-Poincaré method was applied to discretized governing equations to evaluate transverse response of axially moving beams (Pellicano & Zirilli, 1997) and to analyze parametric instability of axially moving elastic beams subjected to multifrequency excitations (Parker & Lin, 2001). The method of normal forms was used to evaluate free vibration of axially moving elastic beams with internal resonances (Pellicano & Zirilli, 1997) as well as forced and parametric vibration of axially moving elastic beams (Pellicano et al., 2001). In (Pellicano & Zirilli, 1997), (Parker & Lin, 2001), and (Pellicano et al., 2001), the mode shapes of axially moving beams were chosen as the base functions and the weight functions, and their orthogonality were employed. The stationary mode shapes can also serve as the base functions and the weight functions to discretize governing equations. Based on the discretization, the Lindstedt-Poincaré method was applied to determine the forced response of axially moving elastic beams (Chen et al., 2007), and the method of multiple scales was applied to evaluate the response of an axially moving viscoelastic beams subjected to multifrequency external excitations (Yang et al., 2009). In their studies, n=2 (Chen et al., 2007) and n=1 (Yang et al., 2009), respectively.

4. Numerical approaches

4.1. Galerkin procedure

Numerical calculation is an effective approach to studying nonlinear vibration of axially moving beams. Based on the numerical solutions of the governing equations, some changing tendencies of vibration characteristics, such as frequencies or amplitudes, with related parameters can be predicted, the approximate analytical results can be verified, and the nonlinear dynamical behaviours can be revealed.

Among numerical approaches, the Galerkin procedure can be applied to discretize the governing equation of nonlinear vibration of axially moving beams. The Galerkin discretization is not only the priority of discretization-perturbations reviewed in Subsection 3.2, but also feasible approach to numerical solutions. Using the 3 order Galerkin discretization of governing equation (in the type of equation (18)) for transverse motion of axially moving viscoelastic beams excited by the changing tension, (Marynowski, 2002) and (Marynowski & Kapitaniak, 2002) numerically investigated the effects of different viscoelastic models, such as the Kelvin model, the Maxwell model, and the standard linear solid model, on the dynamic response and found that different viscoelastic models yield very close numerical results for small damping. The Galerkin procedure has been mainly use to calculate long time nonlinear dynamical behaviors, which will be addressed in Subsection 5.1.

In the application of the Galerkin discretization, the main problem in the actual computations is the complexity of the resulting discretized equations. If the number of terms retained in the Galerkin discretization is rather large, the explicit expression of nonlinear terms is very difficult to obtain. (Chen & Yang, 2006b) proposed a technique to simplify the nonlinear terms in the equations derived from the Galerkin discretization. All nonlinear terms are regrouped to combine the repeated terms and cancel the zero terms. Therefore, the resulting equations can be easily coded for computers and then be effectively calculated. For example, the Galerkin discretization of the governing equation (18) for transverse motion (in the dimensionless form) is

$$\begin{array}{l}{\displaystyle \sum _{i=1}^n{\ddot{q}}_i}\left(t\right)⟨{\phi }_i,{\psi }_k⟩+2\gamma {\displaystyle \sum _{i=1}^n{\dot{q}}_i}\left(t\right)⟨{{\phi }^{\prime }}_i,{\psi }_k⟩+({\gamma }^2-1){\displaystyle \sum _{i=1}^n{q}_i}\left(t\right)⟨{{\phi }^{″}}_i,{\psi }_k⟩\\ +k_f{}^2{\displaystyle \sum _{i=1}^n{q}_i}\left(t\right)⟨{{{\phi }^{″}}^{″}}_i,{\psi }_k⟩+\alpha {\displaystyle \sum _{i=1}^n{\dot{q}}_i}\left(t\right)⟨{{{\phi }^{″}}^{″}}_i,{\psi }_k⟩-\frac{3}{2}{k}_1^2{\displaystyle \sum _{i=1}^n{\displaystyle \sum _{j=1}^n{\displaystyle \sum _{s=1}^n{q}_i\left(t\right)q_j\left(t\right)q_s\left(t\right)⟨{{\phi }^{\prime }}_i{{\phi }^{\prime }}_j{{\phi }^{″}}_s,{\psi }_k⟩}}}\\ -\alpha k_2{\displaystyle \sum _{i=1}^n{\displaystyle \sum _{j=1}^n{\displaystyle \sum _{s=1}^n{q}_i\left(t\right)q_j\left(t\right){\dot{q}}_s\left(t\right)\left(⟨2{{\phi }^{\prime }}_i{{\phi }^{\prime }}_s{{\phi }^{″}}_j,{\psi }_k⟩+⟨{{\phi }^{\prime }}_i{{\phi }^{\prime }}_j{{\phi }^{″}}_k,{\psi }_k⟩\right)}}}=0\left(k=1,2,\mathrm{...},n\right)\end{array}$$E29

If both the base and weight functions are chosen as sine functions, the stationary mode shapes for the simply supported beams, equation (29) can be cast into a from convenient to compute. Evaluating the corresponding inner products, regrouping the nonlinear terms to combine the same terms and canceling all null terms in the resulting equation, one obtains

$$\begin{array}{l}{\ddot{q}}_k\left(t\right)+4{\displaystyle \sum _{j+kisodd}^n\frac{kj}{k^2-j^2}}\left[2\gamma {\dot{q}}_j\left(t\right)+\dot{\gamma }{q}_j\left(t\right)\right]+\left({\gamma }^2-1\right)k^2{\pi }^2{q}_k\left(t\right)+k_f^2{k}^4{\pi }^4{q}_k+\alpha k^4{\pi }^4{\dot{q}}_k\\ =\frac{7k{\pi }^4}{16}\{{\displaystyle \sum _{s=k+1}^{n+k}\left\{\left(s-k\right)q_{l-n}{\displaystyle \sum _{i=\max \left\{1,s-n\right\}}^{\min \left\{s-1,n\right\}}\left[i\left(s-i\right)q_j\left(k_1^2{q}_{s-j}+\alpha k_2{\dot{q}}_{s-j}\right)\right]}\right\}}+{\displaystyle \sum _{l=2}^{n-k}\{\left(n+l\right)q_{k+s}{\displaystyle \sum _{j=1}^{s-1}[j\left(s-j\right)q_j(k_1^2{q}_{s-j}}}\\ +\alpha k_2{\dot{q}}_{s-j})]\}\}-\frac{5k{\pi }^4}{16}{\displaystyle \sum _{s=2}^{k-1}\left\{\left(k-s\right)q_{k-s}{\displaystyle \sum _{i=1}^{s-1}\left[j\left(s-j\right)q_j\left(k_1^2{q}_{s-j}+\alpha k_2{\dot{q}}_{s-j}\right)\right]}\right\}-\frac{k\alpha k_2{\pi }^4}{8}}\{{\displaystyle \sum _{s=k+1}^{n+k}[(s}\\ -k){\dot{q}}_{s-k}{\displaystyle \sum _{i=\max \left\{1,s-n\right\}}^{\min \left\{s-1,n\right\}}j\left(s-j\right)q_j{q}_{s-j}}]+{\displaystyle \sum _{s=2}^{n-k}\left[\left(k+s\right){\dot{q}}_{k+s}{\displaystyle \sum _{i=1}^{s-1}j\left(s-j\right)q_j{q}_{s-j}}\right]}+{\displaystyle \sum _{l=2}^{k-1}\left[\left(k-s\right){\dot{q}}_{n-s}{\displaystyle \sum _{i=1}^{l-1}j\left(s-j\right)q_j{q}_{s-j}}\right]}\}\end{array}$$E30

where the sum is defined to be zero if its lower limit is larger than its upper limit. Although equation (30) seems rather complicated, it is very efficient when used for computer implementing, because almost all repeated nonlinear terms are put together, and terms with zero coefficients are eliminated. In fact, equation (29) contains 2n³ nonlinear terms, while equation (30) contains less than 2n² nonlinear terms. For large n, the difference is significant.

It should be remarked that, based on stationary mode shapes, the even order Galerkin discretization can take the linear gyroscopic terms into full account, while the odd order discretization will miss some effects of the gyroscopic terms.

4.2. Finite difference

The finite difference method is a numerical procedure to solve partial differential equations. The method can be used to discretize both spatial coordinates and time or to discretize spatial coordinates only. In the former case, the procedure consists of four steps: 1 Discretize the continuous spatial domain and temporal interval, on which a partial differential equation is defined, into a discrete finite difference grid; 2 Approximate the individual exact partial derivatives in the partial differential equation by algebraic finite difference approximations; 3 Substitute the finite difference approximations into the partial differential equation in order to derive a set of algebraic finite difference equations; 4 Solve the resulting algebraic finite difference equations.

The finite difference method can be applied to calculated nonlinear vibration of axially moving beams. For example, the method will be employed to solve numerically equation (11) (Chen & Ding, 2010). Introduce the L×T equispaced mesh grid with time step τ and space step h: x_j=jh (j=0, 1,2,…,L, h=l/L); t_n=nτ (n=0,1,2,…,N, τ=T/N), where T is the calculation termination time. Denote the function values u(x,t) and v (x,t) at (x_j,t_n) as uⁿ_j and vⁿ_j. Application of centered difference approximations to the spatial, temporal and mixed partial derivatives leads to

$$\begin{array}{l}u{,}_x=\frac{u_{j+1}^n-u_{j-1}^n}{2h},u{,}_{tt}=\frac{u_j^{n+1}-2u_j^n+u_j^{n-1}}{{\tau }^2},u{,}_{tt}=\frac{u_{j+1}^n-2u_j^n+u_{j-1}^n}{h^2},\\ u{,}_{xt}=\frac{u_{j+1}^{n+1}-u_{j+1}^{n-1}-u_{j-1}^{n+1}+u_{j-1}^{n-1}}{4h\tau }\end{array}$$E31

,and

$$\begin{array}{l}v{,}_x=\frac{v_{j+1}^n-v_{j-1}^n}{2h},v{,}_{xx}=\frac{v_{j+1}^n-2v_j^n+v_{j-1}^n}{h^2},v{,}_{xxx}=\frac{-v_{j-2}^n+2v_{j-1}^n-2v_{j+1}^n+v_{j+2}^n}{2h^3},\\ v{,}_{xxxx}=\frac{v_{j-2}^n-4v_{j-1}^n+6v_j^n-4v_{j+1}^n+v_{j+2}^n}{h^4},v{,}_{xxxxx}=\frac{-v_{j-3}^n+4v_{j-2}^n-5v_{j-1}^n+5v_{j+1}^n-4v_{j+2}^n+v_{j+3}^n}{2h^5},\\ v{,}_{tt}=\frac{v_j^{n+1}-2v_j^n+v_j^{n-1}}{{\tau }^2},v{,}_{xt}=\frac{v_{j+1}^{n+1}-v_{j+1}^{n-1}-v_{j-1}^{n+1}+v_{j-1}^{n-1}}{4h\tau },\\ v{,}_{xxxxt}=\frac{v_{j+2}^{n+1}-4v_{j+1}^{n+1}+6v_j^{n+1}-4v_{j-1}^{n+1}+v_{j-2}^{n+1}-v_{j+2}^{n-1}+4v_{j+1}^{n-1}-6v_j^{n-1}+4v_{j-1}^{n-1}-v_{j-2}^{n-1}}{2\tau h^4}\end{array}$$E32

,

Substitution of equations (31) and (32) into equation (11) leads to a set of algebraic equations with respect to uⁿ_j and vⁿ_j that can be solved as under the boundary conditions (8) and (9) for prescribed parameters and initial conditions. Then the resulting grid values uⁿ_j and vⁿ_j are used in the finite difference schemes as an approximation to the continuous solutions u(x,t) and v(x,t) to equation (11). When the external transverse load is a spatially uniformly distributed periodic force, the amplitude of the beam center displacement changes with the force frequency, which is shown in Fig. 2.

The finite difference method was applied to examine the validity of the two nonlinear transverse models (equations (13) and (14)) and to determine the superiority in the sense of approximating the coupled governing equation (11) of planar vibration. For forced vibration of axially moving viscoelastic beams (Chen & Ding, 2010), the steady-state transverse responses of the beam center calculated from the two transverse models are contrasted with the results based on the coupled equations of planar vibration. Qualitatively, the three models predict the same tendencies with the changing parameters. Quantitatively, there are certain differences. In the view of both the center amplitude and the beam shape, the nonlinear integro-partial-differential equation yields the results closer to those from the governing equation of coupled vibration. The similar result was obtained by the finite difference method for response in free vibrations of axially moving elastic beams (Ding & Chen, 2009a).

Figure 2.

The response amplitude changing with the external excitation frequency.

The finite difference method was used to confirm the analytical results of nonlinear transverse vibration of axially moving beams. For free vibration of axially moving elastic beams, (Pellicano & Zirilli, 1997) compared the beam center displacement changing with time via the Lindstedt- Poincaré method and the normal form method with the numerical solutions via the finite difference method, and found that they are in good agreement. For parametric vibration of axially accelerating viscoelastic beams, (Chen & Ding, 2008) compared the stable steady-state response of the beam center via the method of multiple scales with the numerical solutions via the finite difference method, and demonstrated that they have the same qualitative tendencies changing with the related parameters and are quantitatively with rather high precision.

4.3. Differential quadrature

The differential quadrature method, initiated from the idea of integral quadrature, is an efficient discretization technique to seek accurate numerical solutions using a considerably small number of grid points. The method can be used to discretize both spatial coordinates and time or to discretize spatial coordinates only. In later case, the differential quadrature discretization of a partial-differential equation yields a set of differential-algebraic equations via the following four steps. 1 Discretize the continuous spatial domain, on which a partial differential equation is defined, by grid points; 2 Approximate the individual exact partial derivatives in the partial differential equation by a linear weighted sum of all the functional values at all grid points; 3 Substitute the differential quadrature approximations into the partial differential equation to obtain a set of ordinary-differential-algebraic equations; 4 Solve the resulting ordinary-differential-algebraic equations. The two extensively decisive issues in the applications of the differential quadrature method are to choose grid points and to determine the weighting coefficients for the discretization of a derivative of necessary order.

The differential quadrature method can be applied to calculated nonlinear vibration of axially moving beams. Equation (11) is treated as an example to show the application of the differential quadrature method. Introduce N unequally spaced grid points as (Bert & Malik, 1996) and (Shu, 2001)

$$x_i=\frac{1}{2}\left[1-\cos \frac{\left(i-1\right)\pi }{N-1}\right]\left(i=1,2,\dots ,N\right)$$E33

.

The quadrature rules for the derivatives of a function at the grid points yield

$$\begin{array}{l}u{,}_x\left(x_i,t\right)={\displaystyle \sum _{j=1}^N{A}_{ij}^{\left(1\right)}u\left(x_j,t\right)},u{,}_{xx}\left(x_i,t\right)={\displaystyle \sum _{j=1}^N{A}_{ij}^{\left(2\right)}u\left(x_j,t\right)},\\ v{,}_x\left(x_i,t\right)={\displaystyle \sum _{j=1}^N{A}_{ij}^{\left(1\right)}v\left(x_j,t\right)},v{,}_{xx}\left(x_i,t\right)={\displaystyle \sum _{j=1}^N{A}_{ij}^{\left(2\right)}v\left(x_j,t\right)},\\ v{,}_{xxxx}\left(x_i,t\right)={\displaystyle \sum _{j=1}^N{A}_{ij}^{\left(4\right)}v\left(x_j,t\right)},v{,}_{xxxxx}\left(x_i,t\right)={\displaystyle \sum _{j=1}^N{A}_{ij}^{\left(5\right)}v\left(x_j,t\right)}\end{array}$$E34

,

where the weighting coefficients are the expression

$$A_{ij}^{(1)}=\frac{{\displaystyle \prod _{k=1,k\ne i}^N\left(x_i-x_k\right)}}{\left(x_i-x_j\right){\displaystyle \prod _{k=1,k\ne j}^N\left(x_j-x_k\right)}}\left(i,j=1,2,\cdots ,N;j\ne i\right)$$E35

and the recurrence relationship

$$\begin{array}{l}{A}_{ij}^{(r)}=r\left[A_{ii}^{(r-1)}{A}_{ij}^{(1)}-\frac{A_{ij}^{(r-1)}}{x_i-x_j}\right]\left(r=2,3,4,5;i,j=1,2,\cdots ,N;j\ne i\right),\\ A_{ii}^{(r)}=-{\displaystyle \sum _{k=1,k\ne i}^N{A}_{ik}^{(r)}}\left(r=1,2,3,4,5;i=1,2,\cdots ,N\right)\end{array}$$E36

.

Consider the beam with simply supports at both ends (K=0 in equation (9)). Substitution of equation (34) into equation (11) and modification of the weighting coefficient matrices to implement the boundary conditions (Wang & Bert 1993) lead to the ordinary-differential-algebraic equations

$$\begin{array}{l}{\ddot{u}}_j+2\gamma {\displaystyle \sum _{j=1}^N{A}_{ij}^{\left(1\right)}{\dot{u}}_j}+{\displaystyle \sum _{j=1}^N\left(\dot{\gamma }{A}_{ij}^{\left(1\right)}+{\gamma }^2{A}_{ij}^{\left(2\right)}\right)u_j}+\dot{\gamma }=\\ {\displaystyle \sum _{j=1}^N{A}_{ij}^{\left(1\right)}\left\{\frac{\left(1+{\displaystyle \sum _{j=1}^N{A}_{ij}^{\left(1\right)}{u}_j}\right)}{\rho A\left({\epsilon }_j+1\right)}\left[P_0+A\left(E{\epsilon }_j+\eta {\dot{\epsilon }}_j+\eta \gamma {\displaystyle \sum _{j=1}^N{A}_{ij}^{\left(1\right)}{\epsilon }_j}\right)\right]\right\}}+\frac{f_u\left(x_j,t\right)}{\rho A},\\ {\ddot{v}}_j+{\displaystyle \sum _{j=1}^N\left(2\gamma A_{ij}^{\left(1\right)}+\frac{\eta I}{\rho A}{A}_{ij}^{\left(5\right)}\right){\dot{v}}_j}+{\displaystyle \sum _{j=1}^N\left(\dot{\gamma }{A}_{ij}^{\left(1\right)}+{\gamma }^2{A}_{ij}^{\left(2\right)}+\frac{EI}{\rho A}{A}_{ij}^{\left(4\right)}+\frac{\gamma }{\rho A}{A}_{ij}^{\left(5\right)}\right)v_j}\\ =\frac{1}{\rho A}{\displaystyle \sum _{j=1}^N{A}_{ij}^{\left(1\right)}\left\{\frac{{\displaystyle \sum _{j=1}^N{A}_{ij}^{\left(1\right)}{v}_j}}{{\epsilon }_j+1}\left[P_0+A\left(E{\epsilon }_j+\eta {\dot{\epsilon }}_j+\eta \gamma {\displaystyle \sum _{j=1}^N{A}_{ij}^{\left(1\right)}{\epsilon }_j}\right)\right]\right\}}+\frac{f_v\left(x_j,t\right)}{\rho A},(j=2,3,\mathrm{...},N-1),\\ v_1=v_N=u_1=u_N=0\end{array}$$E37

,

which can be numerically solved via the convenient integration routines.

The differential quadrature method was applied to check the validity and the superiority of the two nonlinear transverse models. For free vibration of axially moving elastic beams (Ding & Chen, 2009a), the transverse responses of the beam center calculated from the two transverse models are contrasted with the results based on the coupled equations of planar vibration. The computational investigation leads to the following conclusions: 1 The differences between the two models are both relatively small for not very large vibration; 2 The model differences increase with the vibration amplitude and the axial speed; 3 The integro-partial-differential equation yields better results.

The differential quadrature method was used to validate the analytical results of nonlinear transverse vibration of axially moving beams. (Chen et al., 2009) developed a differential quadrature scheme to verify the approximate analytical results of stable steady-state response in parametric vibration of axially accelerating viscoelastic beams. Figure 3 shows the comparison, in which the solid and dot lines represent the results of the asymptotic analysis method and the differential quadrature method respectively. The amplitudes from both methods are almost coincided, especially near the exact-resonance (σ=0) and in the first resonance. The differential quadrature method was also applied to confirm the analytical results of the stability regions in linear parametric vibration of axially accelerating beams constituted by the Kelvin model (Chen & Wang, 2009) and the standard linear solid model (Wang & Chen, 2009)

Figure 3.

Comparison of analytical and numerical results.

5. Nonlinear dynamical behaviours

5.1. Galerkin discretization

Axially moving beams undergo periodic vibrations in the aforementioned researches. Nonlinear system may exhibit chaos, steady-state response sensitive to initial conditions thus unpredictable after a certain time and recurrent but either periodic or quasiperiodic hence like a random single with the continuous frequency spectrum. Besides, the dynamical behaviors of nonlinear system may change qualitatively at the critical value of the parametric variation, and the qualitative change is termed as bifurcation.

Many investigations on bifurcation and chaos are based on the Galerkin discretization of various transverse models of axially moving beams. For transverse free vibration of accelerating elastic beams in the supercritical regime, based on 1 order Galerkin discretization, (Ravindra & Zhu, 1998) applied Melnikov’s criterion to find out the parameter condition of occurring chaos and performed numerical simulations to show both period-doubling and intermittent routes to chaos. For transverse harmonically forced vibration of axially moving elastic beams in the supercritical regime, based on 8 order Galerkin discretization, (Pellicano & Vestroni, 2002) observed intricate scenario of chaos, including cascade of bifurcations, blue-sky catastrophes and coexisting chaotic and periodic orbits. Actually, they also considered 12 order Galerkin discretization and found that a few number of degree-of-freedom is sufficient to furnish a good spatial representation and to follow the actual dynamical behaviors. For transverse parametric vibration of axially moving viscoelastic beams excited by the time-dependent tension, based on 4 order Galerkin discretization, (Marynowski, 2004) and (Marynowski & Kapitaniak, 2007) observed the inverse period doubling and inverse Hope bifurcation and occurrences of regular and chaotic motions for beams constituted by the Kelvin model and the standard linear solid model respectively. For transverse parametric vibration of axially accelerating viscoelastic beams, based on 4 order Galerkin discretization, (Chen & Yang, 2006b) constructed numerically the bifurcation diagrams in the case that the axial speed perturbation amplitude, the mean axial speed, or the viscosity coefficient is respectively varied while other parameters are fixed. They also calculated the largest Lyapunov exponent from the discretized governing equation. Numerical results show that, with the increasing speed perturbation amplitude, the increasing mean speed, and the decreasing viscosity coefficient, the equilibrium loses its stability and bifurcates into a periodic motion, and the periodic motion becomes chaotic motion via period doubling bifurcation. In addition, the chaotic motion and the periodic motion exchange alternately for the sufficiently large speed

Figure 4.

Bifurcation versus the dimensionless speed fluctuation amplitude.

Figure 5.

Bifurcation versus the dimensionless viscosity coefficient.

perturbation amplitude and mean speed, and for the sufficiently small viscosity coefficient. Figures 4 and 5 show respectively the bifurcation diagrams versus the dimensionless speed fluctuation amplitude and the dimensionless viscosity coefficient.

5.2. Differential quadrature and time series

The differential quadrature method is an effective numerical technique for initial and boundary problems, and it has much higher precision than the few term Galerkin discretization. However, it has not been applied to calculate nonlinear behaviors of axially moving materials until (Ding & Chen, 2009b). They used the differential quadrature method to investigate bifurcation and chaos of an axially accelerating viscoelastic beam constituted by the Kelvin model. Based on the numerical solutions, analysis of the time series yielded the Lyapunov exponent to identify periodic and chaotic motions. Numerical results show that, with the increasing mean axial speed, the equilibrium loses its stability and bifurcates into a periodic motion, and the periodic motion becomes chaotic motion. The chaotic motion and the periodic motion exchange alternately for the sufficiently large mean axial speed and speed perturbation amplitude. Figures 6 and 7 show the Poincaré map and the largest Lyapunov exponent of periodic and chaotic motions respectively.

Figure 6.

Periodic motion of the beam centre.

Figure 7.

Chaotic motion of the beam centre.

6. Energetics and conserved quantity

6.1. Energetics

It is well known that the total mechanical energy in free vibration of an undamped axially stationary elastic beam with pinned or fixed ends is constant. However, many investigations found that the total mechanical energy associated with free vibration of an axially moving elastic beam is not constant even if the beam travels between two motionless supports. (Barakat, 1968) considered the energetics of an axially moving beam and found that energy flux through the supports can invalidate the linear theories of axially moving beams at sufficiently high transporting speed. (Tabarrok, 1974) showed that the total energy of a traveling beam without tension is periodic in time. (Wickert & Mote, 1989) presented the temporal variation of the total energy related to the local rate of change and calculated the temporal variation of energy associated with modes of moving beams. Considering the case that there were nonconservative forces acting on two boundaries, (Lee & Mote 1997) presented a generalized treatment of energetics of translating beams. (Renshaw et al., 1998) examined the energy of axially moving beams from both Lagrangian and Eulerian views and found that Lagrangian and Eulerian energy functionals are not conserved for axially moving beams. (Zhu & Ni, 2000) investigated energetics of axially moving strings and beams with arbitrarily varying lengths. (Chen & Zu, 2004) proposed energetics of axially moving beams with geometric nonlinearity due to small but finite stretching of the beams. Hence the variation of the total mechanical energy is a fundamental feature of free transverse vibration of axially moving beams. However, all aforementioned investigations on energetics and conserved quantities of axially moving beams have only been limited to transverse vibration, in which longitudinal motion is assumed to be uncoupled and thus neglectable. Actually, the energetics of coupled vibration of axially moving elastic strings (Chen, 2006) can be extended to beams.

Assume that the axially moving beam described at the beginning of Subsection 2.1 is elastic (η=0), is without external excitations (f_u=0, f_v=0), and moves in a constant axially speed (γ=c). Consider the total mechanical energy in a specified spatial domain, the span (0, L). The total mechanical energy consists of the kinetic energy of all material particles and the potential energy resulted from the initial tension, the disturbed tension, and the bending moment caused by the beam deflection due to its motion

$$E={\displaystyle \underset{0}{\overset{L}{\int }}\left\{\frac{\rho A}{2}\left[{\left(c+u{,}_t+cu{,}_x\right)}^2+{\left(v{,}_t+cv{,}_x\right)}^2\right]+\left(P_0+\frac{1}{2}EA\epsilon \right)\epsilon +\frac{1}{2}EIv{,}_{xx}^2\right\}dx}$$E38

Then the time-rate of energy change is

$$$$

$$\frac{dE}{dt}={\displaystyle \underset{0}{\overset{L}{\int }}\frac{\partial }{\partial t}\left\{\frac{\rho A}{2}\left[{\left(c+u{,}_t+cu{,}_x\right)}^2+{\left(v{,}_t+cv{,}_x\right)}^2\right]+\left(P_0+\frac{1}{2}EA\epsilon \right)\epsilon +\frac{1}{2}EIv{,}_{xx}^2\right\}dx}$$E39

Interchanging the order of differentiation and integration and inserting u,_tt and v,_tt solved from equation (11), after some mathematical manipulations, one can express the time-rate of energy change in the boundary values

$$\begin{array}{c}\frac{dE}{dt}={\left[\begin{array}{l}\left(c+u{,}_t+cu{,}_x\right)\frac{\left(P_0+EA\epsilon \right)\left(1+u{,}_x\right)}{\sqrt{{\left(1+u{,}_x\right)}^2+v{,}_x^2}}+\left(v{,}_t+cv{,}_x\right)\frac{\left(P_0+EA\epsilon \right)v{,}_x}{\sqrt{{\left(1+u{,}_x\right)}^2+v{,}_x^2}}\\ +EIv{,}_{xx}\left(v{,}_{xt}+cv{,}_{xx}\right)-EIv{,}_{xxx}\left(v{,}_t+cv{,}_x\right)\end{array}\right]|}_0^l\\ -c{\left\{\frac{1}{2}\rho A\left[{\left(c+u{,}_t+cu{,}_x\right)}^2+{\left(v{,}_t+cv{,}_x\right)}^2\right]+\left(P_0\epsilon +\frac{1}{2}EA{\epsilon }^2\right)+\frac{1}{2}EIv{,}_{xx}{}^2\right\}|}_0^l.\end{array}$$E40

Notice that

$$P_u=\frac{\left(P_0+EA\epsilon \right)\left(1+u{,}_x\right)}{\sqrt{{\left(1+u{,}_x\right)}^2+v{,}_x^2}},P_v=\frac{\left(P_0+EA\epsilon \right)v{,}_x}{\sqrt{{\left(1+u{,}_x\right)}^2+v{,}_x^2}}$$E41

are respectively the longitudinal and transverse components of the tension in the beam, EIv,_xx is the bending moment and EIv,_xxx are the shear, while c+u,_t+cu,_x and v,_t+cv,_x are respectively the absolute velocity in the longitudinal and transverse directions and v,_tx+cv,_xx is the absolute angle velocity. Hence the first term in equation (40) stands for the difference of power of the beam tension, the beam bending moment, and the beam shear acting at two ends. Meanwhile,

$$\widehat{E}=\frac{1}{2}\rho A\left[{\left(c+u{,}_t+cu{,}_x\right)}^2+{\left(v{,}_t+cv{,}_x\right)}^2\right]+\left(P_0\epsilon +\frac{1}{2}EA{\epsilon }^2\right)+\frac{1}{2}EIv{,}_{xx}{}^2$$E42

is the total mechanical energy per unit length. Hence the second term in equation (40) stands for the energy change due to the axial motion of the beam. Physically, equation (40) means the change rate of the energy consisting of two parts: the power of the beam tension, moment and shear applying at two ends and the energy variation resulted from the axial motion.

For a beam with the simple support (K=0 in equation (9)) or the fixed ends (K→∞ in equation (9)), equation (40) leads to, respectively,

$$\frac{dE}{dt}=c{\left[\left(EA-\rho A{c}^2\right)\epsilon \left(1+\frac{\epsilon }{2}\right)+EI\left(\frac{1}{2}v{,}_{xx}-v{,}_{xxx}v{,}_x\right)\right]|}_0^L$$E43

,

$$\frac{dE}{dt}=c{\left[\left(EA-\rho A{c}^2\right)u{,}_x\left(1+\frac{u{,}_x}{2}\right)+EIv{,}_{xx}\right]|}_0^L$$E44

.

For an axially stationary beam, c=0. Equation (40) becomes

$$\frac{dE}{dt}={\left[u{,}_t\frac{\left(P_0+EA\epsilon \right)\left(1+u{,}_x\right)}{\sqrt{{\left(1+u{,}_x\right)}^2+v{,}_x^2}}+v{,}_t\frac{\left(P_0+EA\epsilon \right)v{,}_x}{\sqrt{{\left(1+u{,}_x\right)}^2+v{,}_x^2}}+EIv{,}_{xx}v{,}_{xt}-EIv{,}_{xxx}v{,}_t\right]|}_0^L$$E45

.

If the axially stationary beam is with pinned or fixed ends, equation (45) leads to the conservation of the mechanical energy, which is a well known fact.

6.2. Conserved quantity

Although the total mechanical energy of axially moving beams is generally not constant, there does exist an alternative conserved quantity. (Renshaw et al. 1998) presented both Eulerian and Lagrangian conserved functionals for axially moving beams. (Chen & Zu, 2004) generalized their results to nonlinear free vibration of axially moving beams. They adopted the partial-differential equation (a special case of equation (13)) for axially moving beams undergoing nonlinear transverse vibration. (Chen & Zhao, 2005) also present a conserved functional for a beam modeled by an integro-partial-differential equation derived from the quasi-static assumption (a special case of equation (14)). They applied the conserved functional to verify that the straight equilibrium configuration is stable for beams at low axial speed.

Define the functional

$$I={\displaystyle \underset{0}{\overset{L}{\int }}\left\{\frac{\rho A}{2}\left[\left(u{,}_t{}^2-c^{2}u{,}_x{}^2\right)+\left(v{,}_t{}^2-c^{2}v{,}_x{}^2\right)\right]+\left(P_0+\frac{1}{2}EA\epsilon \right)\epsilon +\frac{1}{2}EIv{,}_{xx}{}^2\right\}dx}$$E46

Evaluation of the temporal differentiation by parts yield

$$\begin{array}{c}\frac{dI}{dt}={\displaystyle \underset{0}{\overset{L}{\int }}\{u{,}_t\left[\rho Au{,}_{tt}+2\rho Acu{,}_{xt}+\left(\rho A{c}^2-EA\right)u{,}_{xx}-\left(EA-P_0\right)\frac{\left(1+u{,}_x\right)v{,}_{x}v{,}_{xx}-v{,}_x{}^{2}u{,}_{xx}}{{\left[{\left(1+u{,}_x\right)}^2+v{,}_x^2\right]}^{3/2}}\right]}\\ +v{,}_t\left[\rho Av{,}_{tt}+2\rho Acv{,}_{xt}+\left(\rho A{c}^2-EA\right)u{,}_{xx}+EIv{,}_{xxxx}-\left(EA-P_0\right)\frac{{\left(1+u{,}_x\right)}^{2}v{,}_{xx}-\left(1+u{,}_x\right)v{,}_{x}u{,}_{xx}}{{\left[{\left(1+u{,}_x\right)}^2+v{,}_x^2\right]}^{3/2}}\right]\}dx\\ -\rho Ac{\left[\left(u{,}_t+cu{,}_x\right)u{,}_t+\left(v{,}_t+cv{,}_x\right)v{,}_t\right]|}_0^L+{\left[\frac{\left(P_0+EA\epsilon \right)\left(1+u{,}_x\right)}{\sqrt{{\left(1+u{,}_x\right)}^2+v{,}_x^2}}u{,}_t+\frac{\left(P_0+EA\epsilon \right)v{,}_x}{\sqrt{{\left(1+u{,}_x\right)}^2+v{,}_x^2}}v{,}_t+EIv{,}_{xx}v{,}_{xt}\right]|}_0^L\end{array}$$E47

.

Substitution of equation (11) with η=0, f_u=0, f_v=0, and γ=c into equation (47) leads to

$$\begin{array}{l}\frac{dI}{dt}=-\rho Ac{\left[\left(u{,}_t+cu{,}_x\right)u{,}_t+\left(v{,}_t+cv{,}_x\right)v{,}_t\right]|}_0^L\\ +{\left[\frac{\left(P_0+EA\epsilon \right)\left(1+u{,}_x\right)}{\sqrt{{\left(1+u{,}_x\right)}^2+v{,}_x^2}}u{,}_t+\frac{\left(P_0+EA\epsilon \right)v{,}_x}{\sqrt{{\left(1+u{,}_x\right)}^2+v{,}_x^2}}v{,}_t+EIv{,}_{xx}v{,}_{xt}\right]|}_0^L\end{array}$$E48

.

At a pinned or fixed end, u,_t=0, v,_t=0, v,_xx=0 or v,_x=0 (hence v,_xt=0). Therefore, equation (48) results in dI/dt=0. There exists functional (46) that is conserved under pinned or fixed boundary conditions for beams moving with a constant axial speed c.

The conserved quantity in a mechanical system is not only mathematically the first integral leading to a reduction in the order of the system, but also reflects the physical essence of the system closely related to the symmetries of the system. Therefore, it is theoretically significant to investigate the conserved quantities. The conserved quantity in a mechanical system can be used to check and develop numerical simulation algorithms. It is also useful for stability analysis and controller design.

7. Concluding remarks

Because an axially moving beam is an effective mechanical model that can be used in diverse engineering fields, many research activities in the area have been witnessed. The chapter summarizes some resent works on modeling, analysis and simulations of nonlinear vibrations of axially moving beams. It will remain to be an active research field. There are many promising topics for future researches, including but surely not limited to the follows: (1) modeling slender structures via sophisticated beam theories such as three-dimensional beams or composite beams, (2) incorporating functionally graded, theromviscoelastic or other advanced materials, (3) accounting for aerodynamic forces and heating and other actions coupled with the vibration, (4) considering complex constraints and coupling such as belts in drive systems, (5) developing analytical approaches especially for coupled vibrations and strongly nonlinear vibrations, (6) investigating convergence, consistency, and stability of numerical procedures, (7) exploring energetics of nonlinear and time-dependent beams under general constraint conditions, (8) understanding complicated dynamical behaviors such as global bifurcations, chaos, patterns, and spatio-temporal chaos.

Acknowledgments

This work was supported by the National Outstanding Young Scientists Foundation of China (No. 10725209), the National Natural Science Foundation of China (No. 90816001), the Specialized Research Fund for the Doctoral Program of Higher Education of China (Grant No. 20093108110005), Shanghai Subject Chief Scientist Project (No. 09XD1401700), Shanghai Leading Talent Program, Shanghai Leading Academic Discipline Project (No. S30106), and the Program for Changjiang Scholars and Innovative Research Team in University (No. IRT0844).