# Nonlinear Vibrations of an Axially Moving Beam with Fractional Viscoelastic Damping

Authors: Zhe Zhang; Hongxiang Yang; Zhongjin Guo; Lili Zhu; Weiyan Liu
Journal: Advances in Civil Engineering
 (2022)

---

## Abstract

The nonlinear vibrations of an axially moving viscoelastic beam under the transverse harmonic excitation are examined. The governing equation of motion of the viscoelastic beam is discretized into a Duffing system with nonlinear fractional derivative using Galerkin’s method. The viscoelasticity of the moving beam is described by the fractional Kelvin–Voigt model based on the Caputo definition. The primary resonance is analytically investigated by the averaging method. With the aid of response curves, a parametric study is conducted to display the influences of the fractional order and the viscosity coefficient on steady-state responses. The validations of this study are given through comparisons between the analytical solutions and numerical ones, where the stability of the solutions is determined by the Routh–Hurwitz criterion. It is found that suppression of undesirable responses can be achieved via changing the viscosity of the system.


---

## Body

## 1. Introduction

The model of axially moving systems that play essential roles is always observed in a wide range of engineering devices, such as power transmission belts, magnetic tapes, paper sheets, chains, pipes conveying fluids, aerial tramways, and fiber textiles. Therefore, research on the dynamic behaviors of such systems has been conducted in the past decades and is still of interest today [1–3]. Especially, the problems of axially moving beams and micro/nano scaled beams axially loaded have widely been tackled in the analysis of some aspects, such as free vibration [4], stability analysis [5, 6], discretization approaches [7], modeling techniques [8], different solution methods [9, 10], and nonlinear dynamics [5–20]. In the architectural design and construction industry, the materials that exhibit excellent damping characteristics have been widely utilized to fabricate structures to enhance their performance from the viewpoint of vibration control. A growing body of research activities has therefore been appearing to explore the nonlinear vibrations of axially moving viscoelastic beams [4–18]. To better understand the damping mechanism, some classical constitutive models such as Kelvin–Voigt [5–19] and three-parameter Zener model [20] are adopted to effectively describe the dynamical responses of the viscoelastic materials. Although the classical models contain combinations of elastic and viscous elements, they do not have sufficient parameters to handle the different shapes of the hysteresis loops reflecting the nature of viscoelastic materials and structures. Consequently, the fractional calculus [21] has been introduced in constitutive relations to obtain a satisfactory solution for the real viscoelastic responses of the materials over a large range of frequency [22, 23]. Although there remain some mathematical issues unsolved, the fractional calculus-based modern viscoelasticity problems are becoming the focus of attention [24–32].This research is devoted to investigating the dynamic behavior of an axially moving viscoelastic beam under a transverse harmonic excitation. It is assumed that the material of the beam obeys the Kelvin–Voigt model based on the fractional Caputo definition. Remarkably speaking, by one-term Garlerkin’s technique, the governing equation of the beam is discretized into a nonlinear Duffing type equation that has the nonlinear fractional operator [26–29]. The first order averaging method is utilized to derive the modulation equations governing the steady state amplitudes and phases of the system. Then, the stability of the solution is studied by the Routh–Hurwitz criterion [33]. Finally, the results of representative calculations are described and briefly discussed from the vibration control point of view.


## 2. Equation of Motion

The present study considers a uniform axially moving viscoelastic beam, shown in Figure1, with density ρ, cross-sectional area A, and moment of inertial I. The initial tension is represented by P0. The beam travels at a constant axial speed v between two motionless ends separated by distance L and is subjected to an external transverse force Fx,T. Here, it is assumed that the excitation is spatially uniform and temporally harmonic: FX,T=FcosWT, where T and X represent the time and the axial coordinate, respectively. Only the bending vibration described by the transverse displacement UX,T is considered here, and Newton’s second law of motion yields(1)ρA∂2U∂T2+2v∂2U∂X∂T+v2∂2U∂X2=∂∂XP0+AσX,T∂U∂X−∂2MX,T∂X2+FcosWT,where σX,T and MX,T are the disturbed axial stress and bending moment, respectively.Figure 1
Schematic of an axially moving beam.The viscoelasticity of the beam material obeys the fractional Kelvin–Voigt model [27, 28, 30] and its constitutive relationship is given as follows:(2)σX,T=E0εLX,T+E1DCαεLX,T.In whichE0 represents the Young’s modulus and E1 denotes the viscoelastic coefficient. To reveal the geometric nonlinearity owing to the small but finite stretching of the beam, we adopt the Lagrangian strain εLX,T, which is defined by(3)εLX,T=12∂UX,T∂X2.In equation (2), DCα denotes the α-order fractional differentiation operator with respect to time T in the Caputo sense given by [21].(4)DCαf=1Γm−α∫0TfmτT−τα+1−mdτ,m−1<α<m,dmdTmf,α=mapositiveinteger,where Γ is the well-known gamma function. Moreover, for a slender beam, the linear moment-curvature relation is adopted.(5)MX,T=E0+E1DCαI∂2UX,T∂X2.Substitution of equations (2), (3) and (5) into equation (1), yields the governing equation of transverse motion of the axially viscoelastic beam(6)ρA∂2U∂T2+2v∂2U∂X∂T+v2∂2U∂X2−P0∂2U∂X2+E0I∂4U∂X4+E1IDCα∂4U∂X4=32AE0∂U∂X2∂2U∂X2+12AE1∂2U∂X2DCα∂U∂X2+AE1∂U∂XDCα∂U∂X⋅∂2U∂X2+FcosWT.In the present investigation, the boundary conditions are taken to be simply supported, i.e.,(7)U0,T=UL,T=0,∂2U0,T∂X2=∂2UL,T∂X2=0.Using the following dimensionless scheme,(8)x=XL,t=TLP0ρA,ω=WLρAP0,u=UL,γ=vρAP0,kf=E0IP0L2,η=E1IP0L2P0ρAL2α,k1=AE0P0,k2=AE02P0P0ρAL2α,p=πFL4P0.Equation (6) becomes(9)∂2u∂t2+2γ∂2u∂x∂t+γ2−1∂2u∂x2+kf2∂4u∂x4+ηDCα∂4u∂x4=32k12∂u∂x2∂2u∂x2+k2∂2u∂x2DCα∂u∂x2+2∂u∂xDCα∂u∂x⋅∂2u∂x2+4πpcosωt.With boundary conditions(10)u0,t=u1,t=0,∂2u0,t∂x2=∂2u1,t∂x2=0.


## 3. Steady State Responses and Stability Analysis

In this study, one term Galerkin’s method is applied to discretize equation (9), and its solution is assumed to be as follows:(11)ux,t=qtsinπx,where qt is the modal coordinate of the beam. Then, a nonlinear fractional ordinary differential equation is derived as follows:(12)q¨+π4ηDCαq+34k2π4qDCαq2+π2π2kf2+1−γ2q+38k12π4q3=pcosωt.For simplicity in the following analysis, equation (12) is rewritten as follows:(13)q¨+μ1DCαq+μ2qDCαq2+ω02q+knq3=pcosωt.With(14)μ1=π4η,μ2=34k2π4,ω0=ππ2kf2+1−γ2,kn=38k12π4.Remarkably, equations (12) or (13) is a nonlinear Duffing type equation having a nonlinear fractional operator DCαq2 [26–28].In what follows, the primary resonance of the fractional system described by (13) will be investigated by the linear averaging method. We mainly focus on examining the effects of the damping parameters and the fractional order on the steady state response of the beam. The stability of the stationary solutions is examined by the Routh–Hurwitz criterion [33].To obtain the primary resonance of the fractional oscillator (equation (13)), we introduce a small parameter ε. We detune the primary resonance by setting ω02=ω2+εσ. Then, equation (13) is rewritten as follows:(15)q¨+ω2q=−εμ1DCαq+μ2qDCαq2+σq+knq3−pcosωt,where the parameter ε merely serves to indicate the assumed smallness of the terms without any physical meaning. The solution at ε=0 can be written in the form as follows:(16)q=Acosθ,θ=ωt+φ,where A and φ are integration constants determined by the initial conditions. When ε≠0 the approximate solution near the primary resonant frequency region can be viewed as a perturbation of the solution given by (16). In this case, A and φ are slowly varying functions of time t such that their derivatives are of oε. According to the averaging method, we assume that the solution of the system takes the following form:(17)q=Acosθandq˙=−Aωsinθ,which implies(18)A˙cosθ−Aφ˙sinθ=0.Differentiating the second equation in equation (17) with respect to t gives the following equation:(19)q¨=−ωA˙sinθ+Aφ˙cosθ−ω2q.Substituting equations (17) and (19) into equation (14) leads to the following equation:(20)A˙=εωsinθGx,t,φ˙=εAωcosθGx,t,where G=μ1DCαq+μ2qDCαq2+σq+knq3−pcosωt.Based on the averaging method, the right-hand side of equation (20) can be replaced by its averages over one period when G is a periodic function with respect to t. However, owing to the presence of the fractional derivative in equation (20), an infinite interval should be considered, and we have the following equation:(21)AverDCαqmcosθsinθ=limT⟶∞1T∫0TDCαqmcosθsinθdθ,m=1,2.Since we consider the steady state only, the fractional order derivative ofq can be easily derived approximately based on the fractional definition given in equation (4) [21, 27, 32], as follows:(22)DCαq≈ωαAcosθ+απ2,qDCαq2≈142ωαA3cosθ+απ2.Subsequently, the average in equation (21) can also be integrated over one period, when the fractional derivative of q is replaced by the respective approximation in equation (22). Then, we have the averaged equation (20) as follows:(23)A˙≈ε2ωπ∫02πsinθG˜A,θdθ,φ˙=ε2ωπA∫02πcosθG˜A,θdθ,where(24)G˜=σAcosθ+Aωαμ1+2α−2μ2A2cosθ+απ2+2α−2μ2ωαA3cos3θ+απ2+kn4A3cos3θ−pcosωt.A direct calculation yields(25)A˙=−ε2ωf11A+f12A3+psinφ,φ˙=ε2ωAf21A+f22A3−pcosφ,where(26)f11=μ1ωαsinαπ2,f12=2α−2μ2ωαsinαπ2,f21=σ+μ1ωαcosαπ2,f22=2α−2μ2ωαcosαπ2+34kn.The small parameterε is reduced in equation (25) by introducing a slow time scale t˜=εt, and we have the following equation:(27)A′=−12ωf11A+f12A3+psinφ,φ′=12ωAf21A+f22A3−pcosφ.where A′=dA/dt˜. The steady state conditions A′=φ′=0 give the following equation:(28)A2f11+f12A22+f21+f22A22=p2,tanφ=−f11+f12A2f21+f22A2−1.Then, the variation of the amplitudeA and phase φ of the steady state primary response as a function of the fractional order. External excitation as well as other control parameters can be determined by these two equations.In what follows, the stability of the solutions of the system in the neighborhood of the equilibrium state is approximately analyzed by exploring the eigenvalues of the Jacobian matrix of equation (27) evaluated at the fixed points of interest. In order to do that, the small disturbances of the response amplitude and phase ΔA=A−A¯ and Δφ=φ−φ¯, are introduced and substituted into equation (27). Then, the linearized equation is obtained as follows:(29)dΔAdtdΔφdt=JA¯,φ¯ΔAΔφ,where(30)JA¯,φ¯=∂F1∂A∂F1∂φ∂F2∂A∂F2∂φA¯,φ¯,with(31)∂F1∂AA¯,φ¯=−12ωf11+3f12A¯2,∂F1∂φA¯,φ¯=−A¯2ωf21+f22A¯2,∂F2∂AA¯,φ¯=12ωA¯f21+3f22A¯2,∂F2∂φA¯,φ¯=−12ωf11+f12A¯2.The eigenfunction of the linearized equation (29) is as follows:(32)λ2+aλ+b=0,with(33)a=−1ωf11+2f12A¯2,b=14ω2f11+3f12A¯2f11+f12A¯2+f21+f22A¯2f21+3f22A¯2.From the Routh–Hurwitz criterion [33], the steady-state response is asymptotically stable if and only if the real parts of the eigenvalues are negative. It can be obtained by the following equation:(34)a<0andb>0.


## 4. Verifications and Numerical Simulations

In the following, a parametric investigation in primary resonance conditions has been conducted to reveal the influences of the control variables on the steady state responses of the beam from the vibration control point of view. The approximate results have been achieved and illustrated by frequency- and forcing amplitude-response curves, in which the solid lines stand for the stable results and the dotted ones correspond to the unstable results. The distinction between the stable solutions and unstable ones is determined by the inequalities (equation (34)) using the Routh–Hurwitz criterion [33]. In calculations, the system parameters in equation (12) have been taken as follows [14]: k1=33.526, k2=0.05, kf=0.173, and γ=0.6.


### 4.1. Effects of Fractional Order on System Response

In Figures2 and 3, the forcing frequency ω near the natural frequency (ω0≈3.038) has been chosen as a bifurcation parameter. The forcing amplitude-response curves are depicted in Figure 4.Figure 2
Forcing frequency-response curves forη=0.04 and p=0.6.Figure 3
Frequency-response curves forp=0.6 and α=0.8.Figure 4
Forcing amplitude-response curves forη=0.02 and ω=4.6.The influence of the fractional orderα on the dynamic behavior of the viscoelastic system is firstly investigated. The response curves for several cases: α=0.6,0.8,and1.0 are plotted with a frequency range near the natural frequency of the system in Figure 2. It is seen apparently that when α=0.6, as the forcing frequency grows gradually from 4 the stable responses increase until a saddle node bifurcation occurs at Pω≈9.592 and another saddle node happens at Qω≈7.547, by which a bistable interval or the hysteresis area can be calculated. For this bistable interval, there are two stable attractors and one unstable attractor in between. This is a typical characteristic of Duffing oscillator due to its cubic nonlinearity in primary resonance [33]. As the forcing frequency increases gradually from point P, the response jumps down from the resonance branch to the nonresonance branch and the system experiences a saddle node bifurcation at point P. Correspondingly, starting at a high forcing frequency on the nonresonance branch, the response undergoes a conversion to the resonance branch again through a jump up at another saddle node bifurcation point Q. In the response curve for α=0.8, although the jump phenomena exist, the positions of the jumps P′ω≈7.192 and Q′ω≈7.096 are shifted to the left, and the width of the jump region (i.e., the hysteresis area) is shrunk obviously. Additionally, the response amplitudes of the system are attenuated at the same time. With the further increase of the fractional order, the hysteresis domain decreases. Eventually, as α increases beyond a certain critical value, the jumps can be eliminated, which is manifested in Figure 2 when α=1. It means that the fractional order can stabilize the system.In addition, the amplitude-response curves are also depicted in Figure4 to verify the above conclusions, where the response amplitudes are plotted as the function of forcing amplitude. Similar conclusions can easily be drawn.


### 4.2. Effects of Viscoelastic Coefficient on the Dynamic Behaviors

The influence of the viscoelastic coefficientη on the dynamic behavior of the viscoelastic system is then investigated and graphically presented by the frequency-response curves for different coefficients, that is α=0.5,0.7,and0.9 in Figure 3. It is found that the influences of the parameter η are much similar to those of the fractional order α. When η=0.02 for example, with the excitation frequency is gradually growing, the amplitude of the stable response continuously rises until the first limit point Pω≈9.821, at which the solution becomes unstable through a saddle node bifurcation. Contrarily, as the forcing frequency decreases, the amplitude exhibits a saddle node bifurcation at the jump-up point Qω≈7.288. Moreover, the hysteresis area is shrunk and eventually eliminated as the viscoelastic coefficient η increases. At the same time, the response amplitudes are attenuated for the cases having a relatively large viscoelastic coefficientη.


### 4.3. Verifications

To verify the previous discussions, the analytical results have been compared with the numerical solutions achieved from the direct numerical integration technique. The primary resonance of the system is shown in Figure5, in which the control parameter-set adopted is chosen as η=0.04, k1=33.526, k2=0.05, kf=0.173, γ=0.6, α=0.6, and p=0.6. Depicted in Figure 5, the stability of the dynamic solutions, such as the “stability limit,” is shown by the conditions that are given in (34). Also, one can observe that the analytical results achieve good agreement with the solutions by the numerical technique.Figure 5
Comparisons of the analytical and numerical solutions whenp=0.6. The solid and dotted lines stand for the stable analytical and the unstable analytical solutions, respectively, and the asterisk marks the numerical result.


## 4.1. Effects of Fractional Order on System Response

In Figures2 and 3, the forcing frequency ω near the natural frequency (ω0≈3.038) has been chosen as a bifurcation parameter. The forcing amplitude-response curves are depicted in Figure 4.Figure 2
Forcing frequency-response curves forη=0.04 and p=0.6.Figure 3
Frequency-response curves forp=0.6 and α=0.8.Figure 4
Forcing amplitude-response curves forη=0.02 and ω=4.6.The influence of the fractional orderα on the dynamic behavior of the viscoelastic system is firstly investigated. The response curves for several cases: α=0.6,0.8,and1.0 are plotted with a frequency range near the natural frequency of the system in Figure 2. It is seen apparently that when α=0.6, as the forcing frequency grows gradually from 4 the stable responses increase until a saddle node bifurcation occurs at Pω≈9.592 and another saddle node happens at Qω≈7.547, by which a bistable interval or the hysteresis area can be calculated. For this bistable interval, there are two stable attractors and one unstable attractor in between. This is a typical characteristic of Duffing oscillator due to its cubic nonlinearity in primary resonance [33]. As the forcing frequency increases gradually from point P, the response jumps down from the resonance branch to the nonresonance branch and the system experiences a saddle node bifurcation at point P. Correspondingly, starting at a high forcing frequency on the nonresonance branch, the response undergoes a conversion to the resonance branch again through a jump up at another saddle node bifurcation point Q. In the response curve for α=0.8, although the jump phenomena exist, the positions of the jumps P′ω≈7.192 and Q′ω≈7.096 are shifted to the left, and the width of the jump region (i.e., the hysteresis area) is shrunk obviously. Additionally, the response amplitudes of the system are attenuated at the same time. With the further increase of the fractional order, the hysteresis domain decreases. Eventually, as α increases beyond a certain critical value, the jumps can be eliminated, which is manifested in Figure 2 when α=1. It means that the fractional order can stabilize the system.In addition, the amplitude-response curves are also depicted in Figure4 to verify the above conclusions, where the response amplitudes are plotted as the function of forcing amplitude. Similar conclusions can easily be drawn.


## 4.2. Effects of Viscoelastic Coefficient on the Dynamic Behaviors

The influence of the viscoelastic coefficientη on the dynamic behavior of the viscoelastic system is then investigated and graphically presented by the frequency-response curves for different coefficients, that is α=0.5,0.7,and0.9 in Figure 3. It is found that the influences of the parameter η are much similar to those of the fractional order α. When η=0.02 for example, with the excitation frequency is gradually growing, the amplitude of the stable response continuously rises until the first limit point Pω≈9.821, at which the solution becomes unstable through a saddle node bifurcation. Contrarily, as the forcing frequency decreases, the amplitude exhibits a saddle node bifurcation at the jump-up point Qω≈7.288. Moreover, the hysteresis area is shrunk and eventually eliminated as the viscoelastic coefficient η increases. At the same time, the response amplitudes are attenuated for the cases having a relatively large viscoelastic coefficientη.


## 4.3. Verifications

To verify the previous discussions, the analytical results have been compared with the numerical solutions achieved from the direct numerical integration technique. The primary resonance of the system is shown in Figure5, in which the control parameter-set adopted is chosen as η=0.04, k1=33.526, k2=0.05, kf=0.173, γ=0.6, α=0.6, and p=0.6. Depicted in Figure 5, the stability of the dynamic solutions, such as the “stability limit,” is shown by the conditions that are given in (34). Also, one can observe that the analytical results achieve good agreement with the solutions by the numerical technique.Figure 5
Comparisons of the analytical and numerical solutions whenp=0.6. The solid and dotted lines stand for the stable analytical and the unstable analytical solutions, respectively, and the asterisk marks the numerical result.


## 5. Conclusions

By the averaging technique, in the present investigations, the primary resonant response of the axially moving viscoelastic beam system has been analytically investigated from the viewpoint of vibration control. The influences of different system parameters, such as fractional order and the viscoelastic coefficient, on the dynamic responses of the viscoelastic beam are graphically illustrated. From the investigation, the following conclusions have been made:(1)
It has shown that the hysteresis or jump phenomenon occurs owing to the existence of multiple solutions on response curves, while the hysteresis area contracts with the increase of the fractional order or viscoelastic coefficient.(2)
A critical value of these two parameters exists, beyond which the hysteresis region can be eliminated. The increase in them reduces the vibration amplitudes.(3)
Both the fractional order and viscoelastic coefficient are powerful factors in stabilizing the system.



---

*Source: 4637716-2022-08-05.xml*