1 Introduction

The impetus comes from the rapid development of process industry and precision machinery. The transverse vibrations of axially moving structures are widely used for over a century and currently attracting considerable attention. Problems about mechanics of axially moving structures are not only theoretically important but also practical.

This is historically true that the vibration of axially moving structures,such as string,beam,membrane,and plate,has been an important theme. Skutch^[1] published the first paper about axially moving strings in German. After half a century,Sack^[2] performed the first English-language investigation also about axially moving strings.

Stability analysis has a long tradition. Mote^[3] first reviewed the dynamic and stability considerations. Varying speeds have been introduced in the model to study the stability of axially moving structures. Miranker^[4] established the equation of transverse vibration for an axially accelerating string for the first time in 1960. From then on,there are many papers about axially accelerating one-dimensional structures. Ulsoy and Mote^[5] investigated transverse vibrations of axially moving plates for the first time in 1982. They analyzed the coupled transverse and torsional vibration of a band saw blade. There are some publications which investigate axially moving elastic plates^[6-22] and viscoelastic plates^[23-29] under different operation conditions.

However,few papers have been given about the dynamic stability of axially accelerating plates. Lee and Ng^[30] used Hamilton's principle and the assumed mode method to investigate the dynamic stability of a moving rectangular plate with axial acceleration and force perturbations. Yang et al.^[31] studied the stability of an axially moving antisymmetric cross-ply composite plate. Tang and Chen^[32] investigated the parametric resonance of axially accelerating viscoelastic plates.

So far,in all available works about axially accelerating membranes and plates and most publications on axially accelerating strings and beams with only exception^[33-39],the tension is presumed to be independent of the longitudinal coordinate. As a result,the values of the tensions are equal at both ends. The outcome contradicts the fact that the system moves with a nonzero acceleration. Because Newton's second law demands that the acceleration results from a nonvanishing composite force. Therefore,the assumption is only an approximation to make the governing equations mathematically easy to be analyzed and cannot be exactly held. To resolve this inconsequence,longitudinally varying tensions are first drawn into axially moving systems on dynamic stability^[33] and combination and principal parametric resonances^[34] of Euler beams,on Timoshenko beams^[35],and on strings^[36].

The rest of this paper is divided as follows. In Section 2,we introduce the linear models of axially moving viscoelastic plates. The governing equations,the initial conditions,and the associated boundary conditions of coupled planar vibrations are derived. In Section 3,the method of multiple scales is used to solve the transverse vibration. In Sections 4 and 5,the stability conditions in the parametric resonances are investigated. In Section 6,the differential quadrature schemes are used to confirm the results of the method of multiple scales. Section 7 ends the paper with the concluding remarks.

2 Modeling

A uniform horizontal viscoelastic plate travels along the $x$-direction with an axial transport time-dependent speed $\gamma $ ($t$),where $t$ is the time. The mass per unit area is $\rho $. Young's modulus is $ E$. Poisson's ratio is $\mu $. The viscoelastic coefficient is $\eta $. The length is $a$,the width is $b$,and the thickness is $h$ in the $x$-,$y$-,or $z$-direction. The displacements of a generic point of the middle surface of the plate are indicated by $u$,$v$,and $w$ in $x$-,$y$-,and $z$-directions,respectively.

The kinetic energy $T$ is

(1)

The variation of the deformation work $\delta W_{\rm d}$ is

(2)

where $z$ is any height measured from the middle surface,$\sigma_{x}$ and $\sigma _{y}$ are the normal stress components,$\tau_{xy}$ is the shear stress component,$\varepsilon_{x}$ and $\varepsilon _{y}$ are the normal strain components,and $\gamma _{xy}$ is the shear strain component. The strain-displacement relationship is

(3)

The plate is constituted by the Kelvin model^[28-29],namely,

(4)

The stress resultants,which are forces per unit length,and the stress moments,which are moments per unit length,are defined as

(5)

$\sigma _{x2}$ is composed of $\sigma _{x1}$ (the normal stress component in the plate with the axial acceleration but without the vibration) and $\sigma _{x}$ (caused by the deflection) of a generic point of the plate at the distance $z$ from the middle surface.

The variation of the virtual work $\delta W_{\rm e}$ done by the excitation per unit area,i.e.,$f_{u}$,$f_{v}$,$f_{w}$ in $x$-,$y$-,respectively,and $z$-directions,and the dissipative force in the $z$-direction can be derived as

(6)

where $ c_{\rm d}$ is the viscous damping coefficient. The plate model is named after a full guide model and a partial derivative model for $k_{1}=1$ and $k_{1}=0$,respectively.

The generalized Hamilton principle is

(7)

With substitution of Eqs. (1) ,(2) ,and (6) into (7) ,the three-dimensional equations of equilibrium of coupled planar vibrations can be described as

(8a)

(8b)

(8c)

The boundary conditions of the plate are

(9a)

(9b)

(9c)

(9d)

(9e)

(9f)

(9g)

(9h)

(9i)

The initial conditions of the plate are

(10a)

(10b)

(10c)

The longitudinal and the tangential displacements are much smaller than the transverse displacement. Therefore,it can be assumed that $u$,$v$,and $w$ are small quantities,$u$ and $v$ are higher-order small quantities than $w$,namely,$u=O(w^{2})$ and $v=O(w^{2})$. As a consequence,the longitudinal change rate of the axial force per unit length will be dependent on the axial acceleration of the axially moving plate^[33]. Therefore,

(11)

where $N_{x0}$ is the initial tension per unit length. The governing equation can be derived as

(12)

Here,the boundary conditions are four edges simply supported for the sake of simplicity,

(13)

It is observed that $w$,$_{yy}$,$w$,$_{yyt}$,$w$,$_{xx}$,$w$,$_{xxt}$,and $w$,$_{xxx}$ must vanish together. We find that the condition (13) can be rewritten as

(14)

If the support is completely the rigid,then $\alpha =0$. Equation (12) can be expressed as

(15)

If the tension variation is ignored,then the tension can be assumed as $N_{x1}=N_{x0}+\alpha \rho \gamma ^{2}$. Equation (2) can be reduced to

(16)

If the support is completely the rigid,and the tension variation is neglected,then Eq. (12) leads to

(17)

Define the dimensionless variables and parameters as follows:

(18)

where $\varepsilon $ is a small dimensionless parameter that is artificially introduced to serve as a bookkeeping device. It accounts for the fact that $\eta $ and $c_{\rm d}$ are small. It will be set as equal to be one in the next analysis. Substitution of Eq. (18) into Eqs. (12) and (14) leads to

(19)

(20)

3 Linear generating system

The axial speed is presumed to be a small simple harmonic variation about the constant mean axial speed $\gamma _{0}$,

(21)

where $\varepsilon \gamma _{1}$ is the amplitude,and $\omega $ is the frequency of the axial speed fluctuation,all in the dimensionless forms. Substituting Eq. (21) into Eqs. (19) and (20) leads to

(22)

(23)

The solution to Eq. (22) is assumed to be

(24)

where $T_{0}=t$ stands for the fast time scale,and $T_{1}=\varepsilon t$ stands for the slow time scale. Substituting Eq. (24) and the relationship

(25)

into Eqs. (22) and (23) ,the equations at each order are separated,

(26a)

(26b)

(27a)

(27b)

The solution to Eq. (26a) is

(28)

where $A_{mn }$ is a complex and undetermined function. c.c. is the complex conjugate of the procee- ding terms. The $mn$th complex mode function $\psi _{mn}$ is^[22]

(29)

where $C_{1n}$ is a complex and undetermined function. $\beta _{jn } (j=1,2,3,$ and 4) and the $mn$th complex frequency $\lambda _{mn}=\delta _{mn}+{\rm i}\omega _{mn}$ can be solved from the characteristic equation and

(30)

where $\omega _{mn}$ is the natural frequency,and $\delta _{mn}$ is the decay rate.

The $mn$th decay rate $\delta _{mn}$ and the $mn$th natural frequency $\omega _{mn}$ can be calculated for $\kappa =0.5$,$\xi =1$,and $\zeta =1$. Figure 1 demonstrates the variation of the first eight dimensionless decay rates and natural frequencies changing with $\gamma _{0}$. Between the mean axial speeds $\gamma _{0}=0$ and $\gamma _{0}=8.997$,the decay rates keep zero while the natural frequencies decrease. In this range,the system is stable. The first natural frequency is equal to be zero,and the first decay rate begins positive at the critical speed $\gamma _{0}=8.998$,and consequently,the divergence happens,and the system is unstable about the zero equilibrium. After $\gamma _{0}=11.196$,the flutter phenomenon occurs,and the system loses stability for the second time.

4 Stability of summation parametric resonance

The summation parametric resonance may occur,while the axial speed variation frequency $\omega $ approaches to the sum of any two natural frequencies.

4.1 Solvability conditions

$\omega $ is expressed as

(31)

where $\sigma $ is a detuning parameter,and $\omega _{kl}$ and $\omega _{k1l1}$ represent the $kl$th and the $k$1$l$1th natural frequencies. In order to investigate summation parametric resonance,the solution to Eq. (26a) can be described by

(32)

where the possible contribution of the $k$2$l$2th mode not involved in the resonance is introduced. Substituting Eqs. (31) and (32) into Eq. (27a) leads to

(33)

where N.S.T. is non-secular terms of the preceding terms.

In order to derive the solvability condition with the nonhomogeneous boundary conditions,it can be assumed that the solution to Eq. (33) is^[41]

(34)

According to Eq. (34) ,$\varphi _{i}(x)$ and $\psi_i(x) (i=kl$,$k$1$l$1,and $k$2$l$2) have different boundary conditions. With substitution of Eq. (34) into Eq. (33) and equalization of coefficients of exp(i$\omega _{kl}T_{0})$,exp(i$\omega _{k1l1}T_{0})$,and exp(i$\omega _{k2l2}T_{0})$ in the resulting equation,the solvability condition should be satisfied with the distributive law of the inner product for complex functions. The complex variable modulation equations for the amplitude and the phase are

(35)

(36)

(37)

where

(38a)

(38b)

(38c)

It can be numerically verified that $d_{h}$ and $f_{h}$ are positive real numbers,while $c_{h}$ is a complex number. Therefore,the solution to Eq. (37) decays to be zero exponentially. Thus,the $k$2$l$2th mode has actually no effect on the dynamic stability.

4.2 Stability conditions

The equations are changed into an autonomous system by introducing the transformation,

(39)

where $p_{h}$ and $q_{h}$ ($h=1$ and 2) are real functions with respect to $T_{1}$. Substituting Eq. (39) into Eqs. (35) and (36) and separating the resulting equation into real and imaginary parts yield

(40)

The characteristic equation of Eq. (40) is expressed as

(41)

The characteristic equation of the Jacobian matrix is

(42)

with

(43)

The Routh-Hurwitz criterion gives

(44)

The calculations lead to

(45)

Equation (44) means that the stability condition is

(46)

where

(47)

4.3 Numerical demonstration

In this paper,unless otherwise noted,the fixed parameters are selected as follows: $\eta =0.000 1$,$c=0.1$,$\kappa =0.5$,$\xi =1$,$\zeta =1$,$k_{1}=1$,and $\gamma _{0}=3$. Then,the first four natural frequencies are $\omega _{11}=17.707 200 25$,$\omega _{12}=48.105 937 76$,$\omega _{21}=47.368 350 24$,and $\omega _{22}=77.445 898 56$.

Figure 2 shows the stability boundaries for the summation parametric resonance in the $\sigma\gamma _{1 }$-plane for different $\eta $. The solid lines denote $\eta =0$,the dashed lines denote $\eta =0.000 1$,and the dash-dot lines denote $\eta =0.000 2$. It can be found from Fig. 2 that $\eta $ has the effect of raising and narrowing the instability zones. The larger $\eta $ leads to the smaller instability range of $\sigma $ for the given $\gamma _{1}$ and the larger instability threshold of $\gamma _{1}$ for the given $\sigma $.

Figure 3 shows the stability boundaries for the summation parametric resonance in the $\sigma \gamma _{1}$-plane for the different viscous damping coefficients $ c_{\rm d}$. The solid lines denote $c_{\rm d}=0$,the dashed lines denote $c_{\rm d}=0.01$,and the dash-dot lines denote $c_{\rm d}=0.1$. It is revealed that the larger $c_{\rm d}$ leads to the smaller instability range of $\sigma $ for the given $\gamma _{1}$ and the larger instability threshold of $\gamma _{1}$ for the given $\sigma $.

Figure 4 shows the stability boundaries for the summation parametric resonance in the $\sigma \gamma _{1 }$-plane for the different models. The solid lines denote the partial derivative model for $k_{1}=0$,and the dashed lines denote the full guide model for $ k_{1}=1$. It can be found that the full guide model has a smaller instability range of $\sigma $ for the given $\gamma _{1}$ and a larger instability threshold of $\gamma _{1}$ for the given $\sigma $ than the partial derivative model.

The stability boundaries based on different models are compared. The stability boundaries in the $\sigma\gamma _{1 }$-plane for the different models are shown in Fig. 5. The dotted line (Model 1) corresponds to Eq. (12) ,the dash-dot line (Model 2) corresponds to Eq. (15) ,the dashed line (Model 3) corresponds to Eq. (16) ,and the solid line (Model 4) corresponds to Eq. (17) . It is revealed that the longitudinally varying tension has a significant effect on the stability boundaries. The tension variation leads to the large instability range of $\sigma $ for the given $\gamma _{1}$ and the small instability threshold of $\gamma _{1}$ for the given $\sigma $. Therefore,ignoring the longitudinally varying tension leads to the unsafe prediction of instability. Especially,when the longitudinally varying tension is taken into account,it is quite clear that the stability boundaries are not sensitive to $\kappa $. When the longitudinally varying tension is considered,the larger $\kappa $ leads to the smaller instability range of $\sigma $ for the given $\gamma _{1}$ and the larger instability threshold of $\gamma _{1}$ for the given $\sigma $. The trend is reversed when the tension variation is neglected.

The comparison between the results with homogeneous boundary conditions (RHBC) and the results with nonhomogeneous boundary conditions (RNHBC) is shown in Fig. 6. The solid lines denote the stability boundaries of the RNHBC. The dashed lines denote the stability boundaries of the RHBC. It shows that the two results have good agreement. The assumption is reasonable in all available papers before.

5 Stability of principal parametric resonance

The principal parametric resonance may occur,while the axial speed variation frequency $\omega $ approaches to the sum of any two natural frequencies.

5.1 Solvability conditions

$\omega $ is expressed as

(48)

The solutions to Eqs. (26a) and (27a) are,respectively,

(49)

With the same steps,we can obtain

(50)

where

(51a)

(51b)

(51c)

It can be numerically verified that $d_{kl}$ is a positive real number,and $c_{kl}$ is a complex number.

5.2 Stability conditions

The next transformation is introduced,i.e.,

(52)

where $p$ and $q$ are real functions. By substituting Eq. (52) into Eq. (50) ,we can obtain

(53)

The characteristic equation of the Jacobian matrix is

(54)

with

(55)

The Routh-Hurwitz criterion gives

(56)

Then,the stability condition is

(57)

5.3 Numerical demonstration

Figure 7 shows the stability boundaries for the first two principal parametric resonances in the $\sigma \gamma _{1}$-plane for different $\eta $. The solid lines denote $\eta =0.000 1$,the dashed lines denote $\eta =0.000 2$,and the dash-dot lines denote $\eta =0.000 3$. The results of both summation and principal parametric resonances yield the same qualitative trends. Compared with those in Fig. 7,the stability boundaries in the higher order resonance are more sensitive to $\eta $.

The stability boundaries in the $\sigma\gamma _{1 }$-plane for $c_{\rm d}$ are shown in Fig. 8. The solid lines denote $c_{\rm d}=0$,the dashed lines denote $c_{\rm d}=0.01$,and the dash-dot lines denote $c_{\rm d}=0.1$. In Fig. 8,it is revealed that the results of both the summation and principal parametric resonances yield the same qualitative trends.

Figure 9 shows the stability boundaries in the $\sigma \gamma _{1}$-plane for different $k_{1}$. The solid lines denote $k_{1}=0$,and the dashed lines denote $k_{1}=1$. It can be found that the results of both the summation and principal parametric resonances yield the same qualitative trends.

Figure 10 shows the stability boundaries for the first two principal and summation parametric resonances. The dash-dot line $P_{1}$ denotes the first principal parametric resonance,the dashed line $P_{2}$ denotes the second principal parametric resonance,and the solid line $S_{12}$ denotes the summation parametric resonance of the first two modes. Overall,the first principal parametric resonance has the largest instability range of $\sigma $ for the given $\gamma _{1}$ and the smallest instability threshold of $\gamma _{1}$ for the given $\sigma $. The second principal parametric resonance has the opposite qualitative trend. However,the summation parametric resonance has the largest instability range of $\sigma $ for the given $\gamma _{1}$ and the smallest instability threshold of $\gamma _{1}$ for the given $\sigma $ in a small region.

The stability boundaries in the $\sigma \gamma _{1}$-plane for different models are shown in Fig. 11. The longitudinal variation of the tension has a significant effect on the stability boundaries. It is noticed that the larger $\kappa $ leads to the larger instability range of $\sigma $ for the given $\gamma _{1}$ and the smaller instability threshold of $\gamma _{1}$ for the given $\sigma $. The tension variation leads to the larger instability range of $\sigma $ for the given $\gamma _{1}$ and the smaller instability threshold of $\gamma _{1}$ for the given $\sigma $.

The effects of the RNHBC and the RHBC are examined in Fig. 12. The solid lines denote the stability boundaries of the plate with the RNHBC. The dashed lines denote the stability boundaries of the plate with the RHBC. The results of the two resonances yield the same qualitative trends.

6 Numerical confirmation via differential quadrature

Consider the domain $x\in $ [0,1] and $y\in $ [0,1],and the numbers of sampling points are $N_{x}$ and $N_{y}$ in the $x$- and $y$-directions,respectively. In this paper,the method of solving the linear generating system and the form of the sampling points in Ref. ^[32] are the same with this paper. Thus,

(58)

Substituting Eq. (28) into Eq. (58) leads to

(59)

Equation (59) can be rewritten into a matrix form,

(60)

where $M$,$G$,and $K$ denote the mass matrix,the gyroscopic matrix,and the stiffness matrix,respectively. Their dimension is ($N_{x}-2) (N_{y}-2) \times (N_{x}-2) (N_{y}-2) $. $\psi _{mn}$ presents the generalized displacement matrix,and the dimension is ($N_{x}-2) (N_{y}-2) \times 1$.

Equation (60) leads to a generalized eigenvalue problem which can be numerically solved. In the following,the sampling points $N_{x}=N_{y}=17$. Figure 13 presents the comparison of the first four complex frequencies for $\kappa =0.5$,$\xi =1$,and $\zeta =1$. The solid lines denote the first four complex frequencies calculated by the analytical method. The black dots denote that they are calculated by the differential quadrature scheme. It can be found that they have almost the same values.

The original equation is

(61)

For $\eta =0.000 1$,$c=0.1$,$\kappa =0.5$,$\xi =1$,$\zeta =1$,$k_{1}=1$,and $\gamma _{0}=3$,Fig. 14 shows the comparison between the analytical and numerical results of the stability boundaries in different resonances in the $\sigma \gamma _{1}$-plane. The solid lines and the dotted lines denote the results of the method of multiple scales and the differential quadrature scheme,respectively. It is noticed that the numerical calculations and the analytical results have reasonable agreement.

7 Conclusions

This paper deals with the axially moving viscoelastic plates subject to the longitudinally varying tensions. The stability is discussed with the approximately analytical method and the numerical investigation. The stability boundaries are analyzed with the method of multiple scales in the different resonances. The results draw the following conclusions:

(i) The first eight complex frequencies vary with $\gamma _{0}$. The system is stable before the critical speed $\gamma _{0}=8.997$. Between the ranges $\gamma _{0}=8.998$ and $\gamma _{0}=11.196$,the system regains the stability. Beyond $\gamma _{0}=11.196$,the flutter phenomenon occurs.

(ii) The introduction of the longitudinally varying tension makes the instability threshold of $\gamma _{1}$ small for the given $\sigma $,and the instability range of $\sigma $ is large for the given $\gamma _{1}$.

(iii) The larger viscoelastic and viscous damping coefficients have the effects of narrowing and raising the instability zones. The viscoelastic coefficient is more sensitive to the stability boundary.

(iv) In the summation parametric resonance,the larger $\kappa $ leads to the smaller instability range of $\sigma $ for the given $\gamma _{1}$ and the larger instability threshold of $\gamma _{1}$ for the given $\sigma $. However,the trend is reversed in the principal parametric resonances if the tension variation is neglected.

(v) The first four complex frequencies calculated by the analytical and numerical methods are almost accordant. The stability boundaries for the different resonances obtained by the two methods have good agreement.