Stability and Dynamics of Viscoelastic Moving Rayleigh Beams with an Asymmetrical Distribution of Material Parameters

Abstract

In this article, vibration of viscoelastic axially functionally graded (AFG) moving Rayleigh and Euler–Bernoulli (EB) beams are investigated and compared, aiming at a performance improvement of translating systems. Additionally, a detailed study is performed to elucidate the influence of various factors, such as the rotary inertia factor and axial gradation of material on the stability borders of the system. The material properties of the beam are distributed linearly or exponentially in the longitudinal direction. The Galerkin procedure and eigenvalue analysis are adopted to acquire the natural frequencies, dynamic configuration, and instability thresholds of the system. Furthermore, an exact analytical expression for the critical velocity of the AFG moving Rayleigh beams is presented. The stability maps and critical velocity contours for various material distributions are examined. In the case of variable density and elastic modulus, it is demonstrated that linear and exponential distributions provide a more stable system, respectively. Furthermore, the results revealed that the decrease of density gradient parameter and the increase of the elastic modulus gradient parameter enhance the natural frequencies and enlarge the instability threshold of the system. Hence, the density and elastic modulus gradients play opposite roles in the dynamic behavior of the system.

1. Introduction

Axially moving beams are broadly explored in various engineering industries, such as telescopic robotic manipulators, band-saw blades, crane hoist cables, high-speed magnetic tapes, and thread lines in the textile industry. For this reason, numerous researchers have devoted considerable attention to the accurate mathematical modeling and dynamic behavior of these applicable structures [1,2,3,4]. Wickret [5] analyzed the nonlinear vibrations of pre-tensioned traveling beams and explained the influence of initial tension and flexural stiffness on the free vibration response. Ghayesh and Amabili [6] numerically solved the forced nonlinear vibrations of axially moving beams and extracted the bifurcation points and steady-state response of the considered system. Chen and Yang [7] inspected the parametric dynamics of axially moving viscoelastic systems with hybrid supports subjected to a fluctuating axial velocity. In the mentioned paper, the natural frequencies, mode shapes, and critical velocity of the system were obtained for different system parameters. Gou and co-workers [8] modeled the coupled thermo-elastic dynamic of axially moving beams. They deduced that the stability borders of these systems are highly dependent on the coupled thermo-elastic factor and beam translation speed. Oz and Pakdemirli [9] considered the vibration mechanism of moving beams with a time-dependent velocity. Their results disclosed that when the amplitude of velocity fluctuations of the system ascends, the stability regions expand. Employing the artificial parameter method, Yang et al. [10] discussed different aspects of the vibrational behavior of moving beams. Additionally, they acquired the closed-form approximated solutions for the natural frequencies of axially moving beams and determined the dynamic stability of these systems systemically. Vibrational performance and stability conditions of accelerating viscoelastic beams surveyed by Chen et al. [11]. They discovered that in the first subharmonic resonance, decreasing the initial stress and increasing the average axial velocity had the same influence on the instability threshold of the beam. Kiani [12] focused on the transverse free vibration instability of nanobeams with axial motion resting on an elastic substrate. He scrutinized the influence of shear and surface energies on the divergence and flutter boundaries. In contrast to the case of Euler–Bernoulli beams with axial motion, the literature on the Rayleigh and Timoshenko moving beams is limited [13,14,15,16,17]. Within this context, Zhou and Chang [18] studied the vibration, stability, and gyroscopic characteristics of rotating Rayleigh beams. They found that by increasing the rotational inertia of the beam, the stability is reduced. Dehrouyeh-Semnani et al. [19] carried out a numerical investigation on the vibration characteristics of axially moving microbeams based on the Timoshenko beam theory. Their numerical results obtained the size-dependent frequencies and critical velocities of the system for clamped-clamped and supported boundary conditions.

Preventing instability and unwanted vibration is a crucial requirement in the engineering design of axially moving continua. Heretofore, engineers have suggested various active and passive vibration control methods to annihilate the excessive vibrations of axially moving structures [20,21]. For instance, Li et al. [22] addressed the active nonlinear dynamic control of axially moving strings using boundary velocity feedback. They showed that the vibrations of the string could be effectively dampened through negative linear velocity feedback at either the upstream or downstream of the system. Zhang et al. [23] attended the passive and adaptive vibration attenuation of axially moving beams attached to nonlinear energy sinks. Their results indicated that the performance of the nonlinear absorbers could be appropriately promoted by adjusting the mass parameter. In all the studies mentioned above, the constituent materials of the moving systems are homogenous and isotropic. During recent years, with the aid of modulating materials properties, engineers improved the mechanical behavior of the moving systems [24,25,26]. In this sense, employing the multiple scales method, the nonlinear forced vibration of axially moving sandwich viscoelastic beams was studied by Li et al. [27]. They obtained the dynamic response and instability boundaries of the considered beam for various values of the core layer thickness and initial tension. Ghayesh [28] concentrated on the dynamic behavior of the moving laminated composite beams and determined the critical velocity, natural frequencies, and mode shapes of the considered system. Lv et al. [29] numerically and analytically inspected the nonlinear lateral dynamics of viscoelastic sandwich beams with varying axial motion and tension. They concluded that applying higher initial tension in the system diminishes the natural frequencies.

To eliminate the structural limitations, researchers introduced novel composite materials known as functionally graded (FG) materials. FG materials are fabricated by smooth and continuous variations of two or more constituent materials in a preferred direction from one surface to another [30,31,32,33]. Compared to the conventional isotropic and laminated materials, FG materials are utilized broadly in the complex engineering applications due to their supreme mechanical characteristics such as lower stress concentration, higher corrosion, and fracture resistance [34,35]. Although FG materials have great potential in the design and development of advanced axially moving components, there are limited studies in this area. In this regard, Piovan and Sampaio [36] carried out the vibration analysis of sliding beams made of FG materials during deploying and retrieving. They supposed that the material characteristics of the system alter according to a power-law function. Furthermore, they evaluated the influence of the acceleration and material constituents on the dynamic response of the system. Sui et al. [37] assessed the lateral dynamics of FG Timoshenko beams with axial motion. They determined the critical velocity and natural frequencies of the considered beam in terms of the length to thickness ratio and FG power index. Kiani [38] presented an analytical study on the transverse and longitudinal vibrations of FG nanobeams within the framework of nonlocal Rayleigh theory. He described the influence of scale parameter and power FG index on the instability threshold of the system. Recently, Yan et al. [39] accomplished the nonlinear vibration of FG beams with axial motion by considering the influence of geometric nonlinearity and axial force. They reported the conditions of occurrence of subharmonic resonance in the system by utilizing the direct multiscale method.

Based on the authors’ knowledge, the previous studies performed on the dynamics of axially moving FG beams considered the gradation of mechanical characteristics of the system along the thickness direction. Despite the importance of metrical gradation through the axial direction, dynamic analysis of moving AFG Rayleigh beams has not been described in the literature. Motivated by this issue, in the current study, the stability enhancement of supported axially moving Rayleigh and EB beams utilizing AFG materials is comprehensively investigated numerically and analytically. Additionally, an in-depth clarification of the effect of various essential parameters such as rotary inertia factor and axial grading of materials on the dynamic features of the axially moving beams, as well as their physical interpretations are presented. It is supposed that the material characteristics of the moving beam are varying exponentially or linearly in the longitudinal direction. The dynamic equation of the considered system is derived from Hamilton’s principle. Then, the reduced-order equation of the system is obtained by implementing Galerkin method, and an eigenvalue problem is carried out. The instability regions of axially moving AFG beam are detected. To confirm the accuracy of the used approach, a comparison study with the available results in the literature is conducted. Time response, natural frequencies, and critical velocity of the considered system are acquired numerically. Additionally, an exact analytical expression is obtained for the critical divergence velocity of AFG moving Rayleigh beams. Finally, the influence of key factors on the vibrational response and instability boundaries of the system are elucidated.

2. Mathematical Modeling

In Figure 1, a viscoelastic AFG moving beam is depicted. It is supposed that the beam with simply supported boundary conditions moves with a constant velocity V and is subjected to axial tension load P. Length, cross-sectional area. as well as moment of inertia of the beam are denoted by L, A, and I, respectively.

The density, ρ(x), and elastic modulus, E(x), of the beam vary along the axial direction linearly or exponentially as follows:

$$\rho \left(x\right)={\rho }_{0}g\left(x\right)$$

(1)

$$E\left(x\right)=E_{0}f\left(x\right)$$

(2)

where:

$$g\left(x\right)=1+\frac{x}{L}\left({\alpha }_{\rho }-1\right)\mathrm{or}g\left(x\right)=e^{\frac{x\ln \left({\alpha }_{\rho }\right)}{L}}$$

(3)

$$f\left(x\right)=1+\frac{x}{L}\left({\alpha }_E-1\right)\mathrm{or}f\left(x\right)=e^{\frac{x\ln \left({\alpha }_E\right)}{L}}$$

(4)

in which α_ρ and α_E represent the density and elastic modulus gradient parameters, respectively. They can be declared as follows:

$${\alpha }_{\rho }=\frac{{\rho }_L}{{\rho }_0}$$

(5)

$${\alpha }_E=\frac{E_L}{E_0}$$

(6)

It is found that at x = 0, E = E_0, and ρ = ρ₀. Furthermore, E = E_L and ρ = ρ_L at x = L. The potential and kinetic energies of the system can be expressed as [5,18,40]:

$$U=\frac{1}{2}{{\displaystyle \int }}_0^L\left(P{\left(w^{\prime }\left(x,t\right)\right)}^2+M{\left(w^{″}\left(x,t\right)\right)}^2\right)\mathrm{d}x$$

(7)

$$\begin{array}{c}T=\frac{1}{2}{{\displaystyle \int }}_0^L\rho \left(x\right)A\left(V^2+{\left(\dot{w}\left(x,t\right)+V{w}^{\prime }\left(x,t\right)\right)}^2\right)\mathrm{d}x\\ +\frac{1}{2}{{\displaystyle \int }}_0^L\rho \left(x\right)I{\left({\dot{w}}^{’}\left(x,t\right)+V{w}^{″}\left(x,t\right)\right)}^2\mathrm{d}x\end{array}$$

(8)

in which prime and overdot indicate spatial and temporal derivatives, respectively. Additionally, w(x, t) presents the lateral displacement of the system, and M is the bending moment and can be calculated by as follows [41]:

$$M=-E\left(x\right)I{w}^{″}\left(x,t\right)+\zeta \left(I{\dot{w}}^{″}\left(x,t\right)+VI{w}^{‴}\left(x,t\right)\right)$$

(9)

where ζ is the viscoelastic coefficient.

According to Hamilton’s principle, one can write [42,43]:

$$\delta {{\displaystyle \int }}_{t_1}^{t_2}\left(T-U\right)\mathrm{d}t=0$$

(10)

By substituting Equations (7)–(9) into Equation (10), the dynamic equation of motion of the considered system can be acquired as:

$$\begin{array}{c}\rho \left(x\right)A(\ddot{w}+2V{\dot{w}}^{\prime }+V^2{w}^{″})+{\rho }^{\prime }\left(x\right)A(V\dot{w}+V^2{w}^{\prime })-\rho \left(x\right)I({\ddot{w}}^{″}+V{\dot{w}}^{‴})\\ -{\rho }^{\prime }\left(x\right)I({\ddot{w}}^{\prime }+V{\dot{w}}^{″})\\ -V\left(\rho \left(x\right)I\right({\dot{w}}^{‴}+V{w}^{⁗})+2{\rho }^{\prime }\left(x\right)I({\dot{w}}^{″}+V{w}^{‴})\\ +{\rho }^{″}\left(x\right)I({\dot{w}}^{\prime }+V{w}^{″}))-P{w}^{″}+E\left(x\right)I{w}^{⁗}+2E^{\prime }\left(x\right)I{w}^{‴}\\ +E^{″}\left(x\right)I{w}^{″}+\zeta I({\dot{w}}^{⁗}+V{w}^{″‴})=0\end{array}$$

(11)

To obtain the dimensionless dynamic equation of the system, the subsequent dimensionless variables are introduced:

$$\xi =\frac{x}{L}, \eta =\frac{w}{L}, \tau =\frac{t}{T}$$

(12)

where:

$$T=\sqrt{\frac{{\rho }_{0}A{L}^2}{P}}$$

(13)

Substituting the dimensionless parameters of Equations (12) and (13) into Equation (11) leads to the subsequent dimensionless dynamic equation:

$$\begin{array}{c}g\left(\xi \right)(\ddot{\eta }+2v{\dot{\eta }}^{\prime }+v^2{\eta }^{″})+\lambda \left(\xi \right)(v\dot{\eta }+v^2{\eta }^{\prime })\\ -\beta \left(g\left(\xi \right)\right({\ddot{\eta }}^{″}+v{\dot{\eta }}^{‴})+\lambda \left(\xi \right)({\ddot{\eta }}^{\prime }+v{\dot{\eta }}^{″})\\ +v\left(g\left(\xi \right)\right({\dot{\eta }}^{‴}+v{\eta }^{⁗})+2\lambda \left(\xi \right)({\dot{\eta }}^{″}+v{\eta }^{‴})+\theta \left(\xi \right)({\dot{\eta }}^{\prime }+v{\eta }^{″}\left)\right))-{\eta }^{″}\\ +v_f^2(f\left(\xi \right){\eta }^{⁗}+2\gamma \left(\xi \right){\eta }^{‴}+\chi \left(\xi \right){\eta }^{″})+\mu ({\dot{\eta }}^{⁗}+v{\eta }^{″‴})=0\end{array}$$

(14)

in which the overdot and prime denote the differentiation with respect to the dimensionless coordinates. Additionally, the dimensionless parameters in Equation (14) are defined as:

$$\begin{array}{c}v=V\sqrt{\frac{{\rho }_{0}A}{P}}, k_f=\sqrt{\frac{E_{0}I}{P{L}^2}}, \beta =\frac{I}{A{L}^2}, f\left(\xi \right)=f\left(x\right), \mu =\frac{\zeta I}{P{L}^3}\sqrt{\frac{P}{{\rho }_{0}A}}\\ \lambda \left(\xi \right)=g^{\prime }\left(x\right)L, \theta \left(\xi \right)=g^{″}\left(x\right)L^2, \gamma \left(\xi \right)=f^{\prime }\left(x\right)L, \chi \left(\xi \right)=f^{″}\left(x\right)L^2\end{array}$$

(15)

It should be mentioned that β and k_f are the rotary inertia factor and the flexural stiffness of the considered system, respectively. It is worth noting that the governing equation of the system reduces to that of AFG moving EB beam by setting β = 0. Furthermore, if the viscoelastic coefficient of the system is supposed to be zero (μ = 0), the dynamic equation degenerates to that of a moving elastic AFG Rayleigh beam.

3. Discretization Technique

With the aim of discretizing the dynamic equation of the system, the lateral displacement of the beam is considered as [44,45,46]:

$$\eta \left(\xi ,\tau \right)={\displaystyle \sum }_{r=1}^n{\phi }_r\left(\xi \right)q_r\left(\tau \right)$$

(16)

in which q_r, n, and φ_r are dimensionless generalized coordinate, the number of basic functions, and acceptable mode shape for the transverse displacement of the system, respectively. The normalized mode shapes of a simply supported beam are given by [47]:

$${\phi }_r\left(\xi \right)=\sqrt{2}\sin \left(r\mathsf{\pi }\xi \right)$$

(17)

By applying the Galerkin discretization procedure, one can obtain the discretized form of the dynamic equation of the system as follows:

$$M\ddot{q}\left(\tau \right)+C\dot{q}\left(\tau \right)+Kq\left(\tau \right)=\mathbf{0}$$

(18)

where q, M, C, and K are the vector of generalized coordinates, mass, damping, and stiffness matrices, respectively, and are given by:

$$q={\left[q_1,q_2,\dots ,q_n\right]}^T$$

(19)

$${\mathbf{M}}_{sr}={{\displaystyle \int }}_0^1\left(g\left(\xi \right){\varphi }_s\left(\xi \right){\varphi }_r\left(\xi \right)-\beta \left(g\left(\xi \right){\varphi }_s\left(\xi \right){\varphi }_r^{″}\left(\xi \right)+\lambda \left(\xi \right){\varphi }_s\left(\xi \right){\varphi }_r^{″}\left(\xi \right)\right)\right)\mathrm{d}\xi$$

(20)

$$\begin{array}{cc}{\mathbf{C}}_{sr}=& v{{\displaystyle \int }}_0^1(2g\left(\xi \right){\varphi }_s\left(\xi \right){\varphi }_r^{\prime }\left(\xi \right)+\lambda \left(\xi \right){\varphi }_s\left(\xi \right){\varphi }_r\left(\xi \right)\\ & -\beta (2g\left(\xi \right){\varphi }_s\left(\xi \right){\varphi }_r^{⁗}\left(\xi \right)+3\lambda \left(\xi \right){\varphi }_s\left(\xi \right){\varphi }_r^{″}\left(\xi \right)+\theta \left(\xi \right){\varphi }_s\left(\xi \right){\varphi }_r^{\prime }\left(\xi \right)))\mathrm{d}\xi \\ & +\mu {{\displaystyle \int }}_0^1{\varphi }_s\left(\xi \right){\varphi }_r^{⁗}\left(\xi \right)\mathrm{d}x\end{array}$$

(21)

$$\begin{array}{cc}{\mathbf{K}}_{sr}=& v^2{{\displaystyle \int }}_0^1(g\left(\xi \right){\varphi }_s\left(\xi \right){\varphi }_r^{″}\left(\xi \right)+\lambda \left(\xi \right){\varphi }_s\left(\xi \right){\varphi }_r^{\prime }\left(\xi \right)\\ & -\beta (g\left(\xi \right){\varphi }_s\left(\xi \right){\varphi }_r^{⁗}\left(\xi \right)+2\lambda \left(\xi \right){\varphi }_s\left(\xi \right){\varphi }_r^{‴}\left(\xi \right)+\theta \left(\xi \right){\varphi }_s\left(\xi \right){\varphi }_r^{″}\left(\xi \right)))\mathrm{d}\xi \\ & +k_f^2{{\displaystyle \int }}_0^1(f\left(\xi \right){\varphi }_s\left(\xi \right){\varphi }_r^{⁗}\left(\xi \right)+2\gamma \left(\xi \right){\varphi }_s\left(\xi \right){\varphi }_r^{‴}\left(\xi \right)\\ & +\chi \left(\xi \right){\varphi }_s\left(\xi \right){\varphi }_r^{″}\left(\xi \right))\mathrm{d}\xi +\mu v{{\displaystyle \int }}_0^1{\varphi }_s\left(\xi \right){\varphi }_r^{″‴}\left(\xi \right)\mathrm{d}\xi -{{\displaystyle \int }}_0^1{\varphi }_s\left(\xi \right){\varphi }_r^{″}\left(\xi \right)\mathrm{d}\xi \end{array}$$

(22)

4. Stability Examination

The second-order differential equations of Equation (18) can be reduced to the following first-order ones as:

$$B\dot{Z}\left(\tau \right)+EZ\left(\tau \right)=\mathbf{0}$$

(23)

where:

$$B=\left[\begin{array}{cc}\mathbf{0}& M\\ M& C\end{array}\right], E=\left[\begin{array}{cc}-M& \mathbf{0}\\ \mathbf{0}& K\end{array}\right], Z\left(\tau \right)=\left[\begin{array}{c}q\left(\tau \right)\\ \dot{q}\left(\tau \right)\end{array}\right]$$

(24)

Assuming Z(τ) = Ae^iωτ yields the following eigenvalue problem:

$$DA-\mathrm{i}\omega I=\mathbf{0}$$

(25)

where I indicates the unity matrix and D = −B⁻¹E. Moreover, ω is the complex-valued natural frequency of the AFG moving Rayleigh beams and can be determined in terms of different key factors such as rotary inertia factor, flexural stiffness, and density and elastic modulus gradient parameters. The stability of the AFG moving Rayleigh beams is profoundly affected by the sign of imaginary and real parts of the natural frequency. Imaginary and real parts of the natural frequency are related to damping and the frequency of oscillation of the system, respectively. When the real part of one of the frequency branches becomes zero while its imaginary part is negative (Re(ω) = 0, Im(ω) < 0), the instability via a pitch-fork bifurcation occurs in the system. The minimum velocity at which the divergence instability happens is recognized as the critical velocity. Furthermore, when the imaginary part of at least one of the frequency branches is negative while its real part is positive (Real(ω) > 0, Image(ω) < 0), the system experiences flutter instability via a Hopf bifurcation [9].

5. Results and Discussion

Firstly, in the case of isotropic system, the results are acquired and compared with the literature. Then, the effect of density and elastic modulus gradations separately and simultaneously on the vibrational frequencies, the dynamic response, and the stability regions of AFG moving Rayleigh and EB simply supported beams are explored and highlighted. It should be mentioned that the dynamic response of the system can be determined numerically by employing the fourth-order Runge–Kutta technique.

5.1. Model Verification

To examine the correctness of the present solution, the real part of the fundamental frequency of an isotropic moving EB simply supported beam under a compressive load versus dimensionless axial velocity is plotted in Figure 2 for various values of flexural stiffness. As revealed in this figure, the results of the current investigation are in good agreement with those presented by Oz and Pakdemirli [9]. According to Figure 2, by increasing k_f, the fundamental frequency of the isotropic system rises.

In Figure 3, the real part of the three vibrational frequencies of isotropic moving simply supported Rayleigh beam under tension load is displayed. The obtained outcomes are in a close agreement with the results obtained by Chang et al. [18]. It should be noted that, to obtain the presented results in Figure 2 and Figure 3, nine terms in Equation (16) are considered. To this end, due to the accuracy of the acquired results, applying nine vibrational bending modes is reasonable for the subsequent computations.

5.2. Effect of Elastic Modulus Variation

In Figure 4a,b, real and imaginary parts of the first two vibrational frequencies of the system versus the axial velocity are displayed, respectively, for the exponential and linear variations of elastic modulus when k_f = 0.5. As it is obvious, the vibrational frequencies of the considered system are purely real when the axial velocity is zero. Afterward, by increasing the velocity, the real part of the vibrational frequencies of the system declines gradually, while their imaginary part is still equal to zero. At the critical axial velocity (v_d), the real part of system frequencies vanished, and the system loses its stability and consequently undergoes the divergence phenomenon. The induced divergence instability in moving structures due to ascending the axial velocity is analogous to that of buckling in classical beams under the compression load [48]. By further increasing the axial velocity, the fundamental frequency of the considered structure becomes purely imaginary, while the second natural frequency declines monotonically. Due to gyroscopic effects in the system, at higher velocities, the beam regains its stability again. In other words, the initiation and termination of divergence instability are correlated to the vanishing of real and imaginary parts of the fundamental vibrational frequency, respectively. Eventually, real parts of the two vibrational frequencies merge into each other via a Paidoussis coupled-mode flutter bifurcation while their imaginary parts divide into two branches with negative and positive values. This phenomenon is related to the flutter instability in the system, and the corresponding velocity is known as the flutter velocity (v_f). In fact, in addition to the velocities lower than critical divergence velocity (v < v_d), a narrow range of axial velocity (between the termination point of the divergence instability and the initiation of flutter instability) exists in which the system is stable at this operational velocity range. It is worth mentioning that the system is no longer stable beyond the critical flutter velocity. As a result, the moving beam experiences a stability evolution of “stable – first mode divergence – stable – coupled-mode flutter”.

According to Figure 4a, the real part of system frequencies ascends by increasing the elastic modulus gradient parameter, especially frequencies of the higher modes. Additionally, increasing α_E leads to ascending the critical divergence and flutter velocities of axially AFG moving EB beams. In other words, it is feasible to hinder the occurrence of undesirable divergence phenomenon by increasing the elastic modulus gradient in moving structures. Since α_E has an increasing role in the stiffness matrix; hence, any increment in α_E leads to a stiffer system and also wider stability regions. In other words, increasing α_E induces the stiffness-hardening effect in the system. Another important feature in Figure 4a,b is that the velocity bandwidth corresponding to the divergence and flutter phenomena in the system (v_d < v < v_f) would be expanded slightly by increasing α_E. Moreover, based on these figures, compared with the exponential variation of elastic modulus, the linear variation leads to a more stable system. As demonstrated in Figure 4b, the damping ratio of the system is higher for α_E > 1 and linear variation of elastic modulus. Accordingly, it is possible to determine the instability thresholds and vibrational behavior of the system by fine-tuning of the elastic modulus gradient parameter.

To better understand the dynamic behavior of the considered structure, the time history of the first generalized coordinate of the system in various velocities are displayed in Figure 5 for k_f = 0.5, β = 0, α_ρ = 1, and α_E = 2. Unit static displacement and zero initial velocity for the first mode ($q_1\left(0\right)=0$ and ${\dot{q}}_1\left(0\right)=0$) are considered for the initial conditions of the system. According to Figure 4a, for v = 1, Real(ω) > 0 and Image(ω) = 0; hence, the system is dynamically stable and generates stable harmonic oscillations. By increasing v, the effective stiffness of the system declines due to the centrifugal force effects. For v = 2.5, the real part of the fundamental frequency of the considered system equals zero, and the beam undergoes instability. In this condition, the dynamic response of the system dramatically grows without oscillation, and the static instability occurs in the system. Increasing the velocity to v = 4 leads to vanishing the imaginary part of frequency, and the beam regains its stability. As the velocity increases, when v = 5, the real part of the frequency grows, while the imaginary part of the vibrational frequency becomes zero. As a result, the oscillation amplitude of the system amplified exponentially by the time. In this case, unlike the divergence instability, the AFG moving beam experiences the flutter instability with growing oscillation. Therefore, the flutter instability is more dangerous than the divergence instability for axially moving structures.

To better understand the stable configuration of the system, stability maps of the system in v_d–k_f and v_d–α_E planes are drawn in Figure 6a,b, respectively. The indicated curves in the stability maps separated the stable and unstable regions, in which, above of each curve, the structure is in the divergence condition. The critical divergence velocity of the system in Figure 4 is consistent with Figure 6. According to Figure 6a, the greater the flexural stiffness, the more stable the system becomes. Therefore, one can conclude that increasing the flexural stiffness parameter has a stabilizing effect on the axially moving systems. Additionally, by increasing α_E, the stability regions of the system expanded. Based on Figure 6a, the linear variation of elastic modulus expands the stability regions of the system more than the exponential one. Moreover, the effect of the elastic modulus gradient is more prominent at higher values of k_f, and the stability borders separate from each other.

According to Figure 6b, compared with the conventional isotropic axially moving beam (α_E = 1), the AFG system is more stable when α_E > 1 and ascending α_E promotes the stability of the system. It is evident that the critical divergence velocity of the system increases by k_f, which can be attributed to the stabilizing effects of the flexural stiffness. By approaching the value of elastic modulus gradient parameter to one (isotropic condition), the stability boundaries of the system are close to each other for linear and exponential distributions, and these boundaries separate from each other by increasing or decreasing of α_E. Additionally, except at α_E = 1, the system is more stable for the linear distribution of elastic modulus in comparison with the exponential one.

5.3. Effect of Density Variation

To study the influence of the density gradient parameter on the system dynamics, the evolution of real and imaginary parts of two vibrational frequencies versus the axial velocity is demonstrated in Figure 7a,b for various values of α_ρ. The influence of α_ρ is more tangible in the vibrational frequencies of higher modes. As displayed, the natural frequencies of the system have a descending trend by increasing the density gradient parameter. Hence, the influence of α_ρ and α_E variations on the vibrational behavior of the system are opposite to each other. Furthermore, the stability region of the system shrinks by ascending the density gradient parameter. Moreover, in contrast to the case of variable elastic modulus scrutinized in the former section, in the case of variable density, the critical divergence and flutter velocities of the exponential distribution is higher than that of the linear distribution.

The density gradient parameter is contributed in the mass, damping, and stiffness matrices, which associated with mass-addition effect, gyroscopic effect, and stiffness-hardening effect, respectively. According to Figure 7a,b, one can conclude that the mass-addition effect is dominant in the system. Another important feature in the frequency diagrams of AFG moving beams is that compared with the case of elastic modulus grading, the system experiences a wider range of vibrational frequencies in the case of density gradation. For this reason, from the perspective of designing, the case of density variation is more effective in avoiding the resonance phenomenon. By scrutinizing Figure 4 and Figure 7, it can be observed that axially grading the materials changes the critical velocities of the system, while it does not alter the order and the type of the system bifurcation series. Accordingly, one can conclude that the quantitative values of natural frequencies and critical velocities are strictly dependent on axial grading of materials, but the qualitative stability of the system does not vary by the axial gradation of materials. Based on Figure 4 and Figure 7, it can be mentioned that the vibrational behavior of the axially moving systems is highly dependent on the α_ρ and α_E, as well as the type of their distributions.

To better describe the dynamic behavior of AFG moving beams, the stability maps in v_d–α_ρ and v_d–β planes are plotted in Figure 8a,b, respectively. As demonstrated in Figure 8a, the AFG moving beam is more stable for α_ρ < 1 in comparison with the isotropic one. It is crucial to mention that the stability borders of exponential and linear distributions approach to each other for the density gradient parameters close to one (i.e., the isotropic case). Moreover, the critical velocity of the structure is higher for the exponential distribution of density in comparison with the linear one, especially for the higher and lower values of α_ρ. In other words, the difference between the stability borders of linear and exponential distributions can be magnified by selecting a higher or lower density gradient.

According to Figure 8a,b, since increasing the rotary inertia factor displaces the stability boundaries toward smaller velocities, EB theory predicts less stability for the axially moving structures in comparison with Rayleigh beam theory. Additionally, for any constant values of β, the stability regions of the system shrink by increasing α_ρ. Generally, compared with the stability maps in v_d–k_f and v_d–α_E planes analyzed in the former section, it can be declared that the stability maps in v_d–β and v_d–α_ρ are overall descending with increasing β and α_ρ. This implies that increasing β and α_ρ can destabilize the system and leads to the decrement of the critical divergence velocity of the structure. It should be mentioned that the indicated stability borders in Figure 7 and Figure 8 are in agreement. Based on Figure 7 and Figure 8, the critical velocity of the structure is substantially dependent on α_ρ. This dependency is more pronounced for the higher and lower values of α_ρ.

5.4. Effect of Simultaneous Elastic Modulus and Density Variations

According to previous sections, it can be concluded that the variations of elastic modulus and density along the axial direction of the structure have remarkable effects on the vibrations of axially moving beams. Additionally, it is demonstrated that mitigating the excessive oscillations is feasible by adjusting α_E and α_ρ, separately. Based on Figure 4, Figure 5, Figure 6, Figure 7 and Figure 8, variations of elastic modulus and density along the longitudinal direction of the beam have opposite effects on the stability of the system. As a result, these parameters can provide additional degrees of freedom to modify the dynamic characteristics of axially moving structures. In other words, it is possible to significantly improve the performance of axially moving systems by simultaneous fine-tuning α_E and α_ρ. Therefore, determining the role of simultaneous gradation of the material properties on the stability of the moving structures is of great significance. In this section, stability characteristics of the system are explored by considering the coupled variation of density and elastic modulus through the axial direction (simultaneous mass-addition and stiffness-hardening effects). Additionally, the exact analytical expression is presented for critical divergence velocity. Furthermore, a comparison is carried out between different solution procedures.

When the system experiences the divergence phenomenon, the effective stiffness of the system equals zero. In other words, at the divergence threshold, the lowest vibrational frequency of the structure becomes zero. Consequently, the critical velocity of the system is related to the first mode and obtained by reducing Equation (16) to the following equation:

$$M_{11}{\ddot{q}}_1\left(\tau \right)+C_{11}{\dot{q}}_1\left(\tau \right)+K_{11}{q}_1\left(\tau \right)=0$$

(26)

in which the subscript one represents the first natural mode. Considering linear variations for material properties of the system leads to the following equation:

$$K_{11}={\pi }^2\left(\beta \left({\alpha }_{\rho }+1\right)v_d^2-\left({\alpha }_E+1\right)k_f^2\right)+v_d^2\left({\alpha }_{\rho }+1\right)-2$$

(27)

Equation (27) confirms that the critical velocity of the AFG moving system is dependent on the β, α_E, α_ρ, as well as k_f of the system. To better inspect the stability of the two-dimensional contour plot for the critical divergence velocity in α_E–α_ρ and k_f–β planes are drawn in Figure 9a,b, respectively. As shown in Figure 9a, the critical divergence velocity of the system enhances by increasing the α_E and decreasing α_ρ and vice versa. As illustrated in Figure 9b, the critical divergence velocity of the system increases by decreasing the rotatory inertia factor and increasing the dimensionless flexural stiffness and vice versa. Additionally, it can be deduced that the influence of the flexural stiffness and the rotary inertia factor on the stability of the system is opposite to each other. In other words, contrary to the effect of α_ρ and β, increasing the elastic modulus gradient parameter and dimensionless flexural stiffness enhance the buckling strength of the system. Accordingly, simultaneously selecting larger values of α_E and k_f as well as choosing smaller values of β and α_ρ leads to a more stable system and, consequently, improves the performance of the axially moving structures.

In Figure 10a–d), the instability thresholds of the system when the elastic modulus and density gradient parameters are equal (α_E = α_ρ = α) are demonstrated, and the validity of the analytical expressions are examined. According to these figures, the critical divergence velocities acquired by the numerical approaches are in a close agreement with those acquired analytically. As seen in Figure 10a, increasing the material gradient parameter, α, leads to a slight decrease in the critical divergence velocity of the structure. Therefore, compared with the stability maps in previous sections (Figure 6 and Figure 8), it can be inferred that the density gradation (mass-addition effect) shows a dominant role in the stability of the system, and elastic modulus gradation (stiffness-hardening effect) has less effect on the vibrational behavior of the structure. It should be mentioned that compared with previous sections, as density and elastic modulus vary simultaneously, stability borders are closer to each other and are less sensitive to the axial material gradation when the effect of flexural stiffness is highlighted. According to Figure 10b, the stability regions of the system shrink by increasing the gradient parameter, and the curves have a descending trend. Moreover, in the higher values of dimensionless flexural stiffness, variation in the gradient parameter has a minor influence on the instability threshold of AFG moving beam; this point can be helpful at the design stage. As depicted in Figure 10c,d, the stability regions shrink by increasing the gradient parameter and rotary inertia factor. Compared with the isotropic case, when α < 1, the axially moving beam is more stable. Furthermore, at sufficiently high gradient parameter, the divergence velocity of the structure experiences small changes.

The critical velocity of the system in the three-dimensional space of α, β, and k_f is plotted in Figure 11. It is detected that the instability threshold of moving beams is entirely dependent on the axial material gradation, dimensionless flexural stiffness, and rotary inertia factor. As it is obvious, the variation of the material gradient parameter and the rotary inertia factor have the same effects on the stability boundaries, and by increasing each of these parameters, the divergence strength of the system diminishes. While ascending the flexural stiffness enhances the stability of the structure. Consequently, in the simultaneous presence of higher values of dimensionless flexural stiffness, lower values of the material gradient, and rotary inertia factor, the performance of axially moving systems would be improved. Therefore, these parameters could be introduced as the key factors for the vibration control of moving continua. Generally, AFG moving beams are more flexible to adjust their vibrational behavior in comparison with isotropic ones. As a result, it can be claimed that, compared with isotropic materials, AFG ones have a better performance in axially moving structures.

5.5. Effect of Viscoelastic Material

Finally, to explore the influence of viscoelastic materials on the dynamics of the moving beam, real and imaginary parts of the two vibrational frequencies of the system are plotted versus the velocity in Figure 12a,b, respectively, for μ = 0.001. The critical divergence velocity of the structure does not change by ascending the viscosity coefficient. This point can be verified by the analytical solution presented in the previous section. Since the viscoelastic structure is non-conservative, the vibrational frequencies are complex before the buckling, especially the frequency of higher modes. According to Figure 12b, imaginary parts of the frequency curves lose their symmetry toward the x-axis when the system is viscoelastic. Additionally, it can be understood that the viscoelastic structure experiences the stability evolution of “stable-first mode flutter-second mode divergence”. As a result, it can be established that compared with the isotropic and axially graded systems, utilizing viscoelastic materials change the stability evolution of the axially moving continua. Generally, it can be stated that the qualitative stability of the beams with axial motion is dependent on the effect of viscoelastic materials, while, axial gradation of materials plays a significant role in determining quantitative values of the critical velocity and the natural frequencies of the structure.

6. Conclusions

In this study, structural dynamics and possible vibrational instabilities of AFG moving beams are investigated analytically and numerically in detail. Linear and exponential distributions for the axial grading of material properties are considered for the system. By applying the Galerkin discretizing technique and eigenvalue analysis, the natural frequencies, dynamic response, divergence and flutter instability boundaries of the system are obtained for coupled effects of beam velocity, dimensionless flexural stiffness, viscosity coefficient, and axial material gradation. Closed-form mathematical expression is derived for the critical divergence velocity of the system. Stability maps and two-dimensional contour plots of critical divergence velocity are plotted for Rayleigh and EB beams. The main outcomes of the current study may be narrated as follows:

• Increasing the density/elastic modulus gradient parameter has a destabilizing/stabilizing effect on axially moving beams. Compared with isotropic axially moving beams, the system is more stable when density/elastic modulus decreases/increases along the axial direction.
• In the case of density/elastic modulus variation, exponential/linear distribution leads to a more stable system.
• In the case of simultaneous axial variation of elastic modulus and density, the effect of density gradation on the vibrational configuration of the system is dominant.
• The higher flexural stiffness, and the lower rotary inertia factor, the more stable the structure becomes. Moreover, the influence of axial material gradation on the stability boundaries of the system is more tangible at higher and lower values of flexural stiffness and rotary inertia factor.
• Compared with isotropic and moving axially graded beams, utilizing the viscoelastic material changes the stability evolution of the system.

It is demonstrated that compared with isotropic moving beams, axially moving beams present a better performance by designating appropriate viscoelastic AFG materials.