# Analysis of Complex Modal Characteristics of Fractional Derivative Viscoelastic Rotating Beams

Authors: Tianle Lu; Zhongmin Wang; Dongdong Liu
Journal: Shock and Vibration
 (2019)

---

## Abstract

For the transverse vibration problem of a fractional derivative viscoelastic rotating beam, the differential equation of the system is obtained based on the Euler–Bernoulli beam theory and Hamilton principle. Then, introducing dimensionless quantities to differential equations and boundary conditions, the generalized complex eigenvalue equations of the system are obtained by the differential quadrature method. The effects of the slenderness ratio, the viscoelastic ratio, the hub radius-beam length ratio, and dimensionless hub speed and fractional order on the vibration characteristics of fractional derivative viscoelastic rotating beams are discussed by numerical examples. Numerical calculations show that when the dimensionless hub speed is constant, the real part of complex frequency increases with the increase of the fractional order, and the higher-order growth trend is more obvious. Through the study of displacement response at different points on the beam, it can be seen that the closer to the free end, the larger the response amplitude. And, the amplitude of response has been attenuated, which is also consistent with the vibration law of free vibration considering damping.


---

## Body

## 1. Introduction

With the rapid development of manufacturing industry and mechanical structure engineering, rotating structure is more and more widely used in mechanical engineering, such as internal combustion engines, satellites, turbines, and CNC machine tools.Many scholars have made a lot of research on the dynamics of rotating beams. Bannerjee and Su [1] analyzed the free vibration characteristics of Timoshenko beams by the dynamic stiffness method and compared the results with Euler–Bernoulli theory. In consideration of Coriolis force and centrifugal force, Tian et al. [2] investigated the influence of hub radius, slenderness ratio, and angular velocity on the vibration characteristics of rotating beam. Aksencer and Aydogdu [3] studied the vibration problem of rotating composite beams by using the Ritz method, and the effects of several geometric factors on the vibration characteristics of composite beams under different boundary conditions were analyzed. Considering geometric nonlinearity, Zhang et al. [4] studied the nonlinear vibration of a variable cross-section viscoelastic rotating beam, which is based on Kelvin–Voigt viscoelastic constitutive relation, and the amplitude-frequency response of the rotating beam is obtained by the multiscale method. Li et al. [5] studied the vibration characteristics of a rotating FGM circular section flexible beam. Srivatsa and Ranjan [6] solved the vibration problem of a rotating Rayleigh beam with variable cross-section and obtained the frequency and mode shapes of the beam by the finite element method. Tang et al. [7] investigated the free vibration of the rotating cone cantilever beam based on the Hamiltonian principle, and then analyzed the effects of angular velocity and the moment of inertia on the free vibration of the system. Sofiyev and Hui [8] used FOSDT and Donnell shell theory to study the vibration and stability of cylindrical shells with functionally graded materials under mixed boundary conditions and analyzed the influence of factors such as FG distribution on the critical parameters of FGMCSS containing MBCS. Based on FOSDT and isotropic constitutive relation, Yao et al. [9] studied the nonlinear dynamic response of a rotating cylindrical shell with pretorsion angle by using Hamilton principle, and the effects of excitation parameters, damping coefficient, rotational speed, prerotation angle, and pretorsion angle on the nonlinear dynamic response of the system were analyzed. Based on the dynamic stiffness method, Banerjee and Jackson [10] investigated the free vibration characteristics of the rotating cone-shaped Rayleigh beam and analyzed the effects of the viscosity-elasticity rate on the free vibration characteristics of the rotating conical Rayleigh beam. Based on Euler–Bernoulli, Timoshenko, Reid, and Levison’s four different beam theories, Chakraverty and Behera [11] used the DQM to investigate the buckling of nanobeams and analyzed the effects of nonlocal parameters, aspect ratio, boundary conditions, and nonuniform parameters on the critical buckling load parameters.The fractional derivative viscoelastic model is considered to be a model that can accurately describe the constitutive relation of viscoelastic materials, and it can describe the mechanical behavior of materials in a wide frequency range and determine the experimental parameters required for the model. Based on fractional derivative theory, Gao et al. [12] investigated the viscoelastic mechanical parameters of three standard linear solid material properties. Tarasov [13] investigated the characteristics of fractional derivatives of noninteger order and proposed the nonlocality principle of fractional derivatives. Achouri et al. [14] investigated the initial-boundary value problem of Euler–Bernoulli beam equation with fractional derivative boundary conditions and analyzed the existence and attenuation of its initial-boundary value solution. Based on Euler–Bernoulli beam theory, Ansari et al. [15] studied the geometric nonlinearity vibration of fractional order viscoelastic nanobeams. Bahraini et al. [16] studied the free vibration of viscoelastic cantilever beams with different boundary conditions under concentrated loads and analyzed the influence of fractional derivatives on the large deflection of viscoelastic cantilever beams.At present, the theory and application of the differential quadrature method (DQM) are quite mature, and it is mainly used to solve some problems such as high-order differential equations with variable coefficients. Li and Ying [17] solved the transverse vibration of elastic Euler–Bernoulli beams under interharmonic loads by the differential quadrature method and analyzed its vibration characteristics. Based on the Euler–Bernoulli beam model, Zhang et al. [18] studied the large deformation mechanical behavior of members by the differential quadrature method. Cao [19] investigated the static response of reinforced concrete (RC) beams using the generalized differential quadrature method by virtual work principle. Through the differential quadrature method, Matbuly et al. [20] solved the free vibration problem of functionally graded elastic beams and studied the influence of geometric parameters on the natural frequency of the beam. Li et al. [21] studied the dynamic properties of viscoelastic Timoshenko beams by the differential quadrature method and the nonlinear numerical method. Considering the geometric parameters, moment of inertia, and quadratic coupling deformations of variable cross-section beams, Wang and Li [22] studied the transverse vibration characteristics of a rotating cone cantilever beam by using Hamiltonian variational principle, and the DQM was used to solve the differential equation of motion. Based on the thin plate theory, Zhou and Wang [23] solved the vibration problem of variable thickness viscoelastic plates by the differential quadrature method and studied the influence of length to diameter ratio and thickness to diameter ratio on the vibration of linear thickness viscoelastic plates.In this paper, based on Euler–Bernoulli beam theory and Hamilton principle, the vibration model and motion differential equation of the fractional derivative viscoelastic rotating beam are established. The dimensionless quantity is introduced, and the differential equations and boundary conditions are processed in a dimensionless manner. Then, the differential quadrature method is used to discretize the differential equation and boundary conditions, and the generalized complex eigenequation of the system is obtained. The effects of slenderness ratio, the hub radius-beam length ratio, the viscoelastic ratio, dimensionless hub speed, and fractional order on the complex frequency characteristics of the system are analyzed.


## 2. One-Dimensional Fractional Kelvin–Voigt Model Theory

The classical Kelvin–Voigt model is composed of a spring and a Newton viscous pot in parallel, and the stress-strain relationship of the fractional Kelvin model is as follows:(1)σxt=Eεxt+ηDαεxt,where E and η are the Young’s modulus of elasticity and damping modulus of the material, εxt is the Lagrange strain of finite deformation; Dα is the Riemann-Liouville fractional differential operator, and α is a fractional order.Set0≤α≤1, the α order Riemann–Liouville fractional derivative of the function ft defined on 0,t is expressed as follows [12]:(2)Dαft=f0Γ1−αt−α+1Γ1−α∫0tf′τt−ταdτ.


## 3. Differential Equation of Motion

### 3.1. Deformation Field Description

Figure1(a) shows that the schematic diagram of the fractional derivative viscoelastic rotating hollow beam, which has length L of the beam, radius a, and rotating angular velocity ω of the hub. The oxy and OXY are the floating coordinate system and the inertial coordinate system, respectively. The θ represents rotation angle of the viscoelastic beam in a wide range of motion, that is, the rotation angle of the oxy is relative to the OXY. The ox-axis coincides with the axis of the viscoelastic hollow circular section beam before deformed.Figure 1
Dynamic model of fractional derivative viscoelastic rotating hollow beam. (a) Fractional derivative viscoelastic rotating hollow beam. (b) Viscoelastic rotating hollow beam profile.
(a)
(b)In Figure1(b), the cross-section of the beam is a hollow section with wall thickness h and Rn is the radius of the wall thickness for any cross-section at the middle line. In the n−s coordinate system, n is the normal direction and s is the tangential direction.In Figure2, ra is a radius vector of the hub, xp is the vector diameter of any point P on the axis of a fractional derivative viscoelastic beam relative to the point o in the floating coordinate system when it is not deformed, and up is a displacement vector of the point P.Figure 2
Displacement field at any point on the beam axis.The position vector of pointP relative to the point O can be written as follows:(3)rp=ra+Aθxp+up,where(4)ra=acosθasinθT,xp=x0T,Aθ=cosθ−sinθsinθcosθ,up=upxupy=u+uc+ubw=u−12∫0xw,x2dx−yw,xw,where Aθ, u, and w are the transformation matrices of the oxy relative to the OXY, the axial displacement, and the lateral bending deformation displacement on the beam axis, respectively; ub=−yw,x is the axial displacement caused by the rotation of the beam cross-section, uc=−1/2∫0xw,x2dx is the secondary deformation coupling of the beam, that is, the axial expansion amount caused by the lateral bending deformation of the fractional derivative viscoelastic rotating hollow beam, and “,x” represents the first derivative of x.Taking the first derivative of time for equation (3), we get the velocity vector of point P as follows:(5)r˙p=r˙a+θ˙×Αθxp+up+Αθu˙p=−a+x+uc+ubθ˙sinθ−wθ˙cosθ+u˙c+u˙bcosθ−w˙sinθa+x+uc+ubθ˙cosθ−wθ˙sinθ+u˙c+u˙bsinθ+w˙cosθ,where θ˙=0−θ˙θ˙0 is the antisymmetric matrix of angular velocity θ˙ and “⋅” represents the first derivative of t.


### 3.2. Differential Equations of Motion

The total kinetic energy of the system can be seen as the sum of the kinetic energy of the hub and the kinetic energy of the fractional derivative viscoelastic rotating beam, namely(6)Ts=TH+T=12JHθ˙2+12ρ∫Vr˙pTr˙pdV=12JHθ˙2+13ρIθ˙2∫0Lw,x2dx+13ρI∫0Lw˙,x2dx+12∫0LρAw˙2+u˙c2+θ˙2w2+a+x+uc2dx+∫0LρAθ˙−u˙cw+w˙a+x+ucdx,where JH is the inertia moment of the hub with respect to the central axis.The strain energy of the system can be divided into two parts, one part is potential energy and the other part is dissipative energy.The potential energy can be expressed as follows:(7)V0=∫V12σxεxdV=∫V12Eεx2dV,where σx, εx, and V are the normal stress, normal strains of the beam, and the volume of the viscoelastic beam, respectively.In equation (4), not considering the axial displacement and ignoring nonlinear term, the normal strain can be obtained by the relationship between strain and displacement, namely,(8)εx=upx,x=−12∫0xw,x2dx−yw,x,x=−ys−nz,sw,xx−12w,x2.According to Figure1(b), the geometric relations can be obtained:(9)dzds=−sinϕ=−sinsRn,dyds=cosϕ=cossRn.Then, the potential energy can be expressed as follows:(10)Vo=∫0L∫A12Eεx2tdAdx=∫0L∫−h/2h/2∫0πRn12ERn+nsinsRnw,xx+12w,x22Rn+nRndsdndx=12EI∫0Lw,xx2dx+18EA∫0Lw,x4dx,where I and A are the inertia moment and the cross-sectional area of the beam, respectively,(11)I=164π2Rn+h4−2Rn−h4=14π4Rn3h+Rnh3,A=πRn+h22−Rn−h22=2πRnh.The expression of the system’s dissipative energy is(12)Wc=∫V∫0εxtηDαεxtdεxtdV.From equation (2), the following equation can be obtained(13)Dαεxt=εx0Γ1−αt−α+1Γ1−α∫0tε˙xτt−ταdτ.Substituting equation (12) into equation (13), the following equation can be obtained:(14)Wc=∫VηDαεxtεxtdV=∫V∫0εxtηεx0Γ1−αt−αdεxtdV+∫V∫0εxtη1Γ1−α∫0tε˙xτt−ταdτdεxxdV=ηI∫0Lw,xxDαw,xxdx+14ηA∫0Lw,x2Dαw,x2dx.In this paper, based on Hamiltonian variational principle, the differential equation of fractional derivative viscoelastic rotating beam is derived. Hamilton’s variational principle takes the form of(15)δ∫t1t2Ldt=δ∫t1t2Ts−V0−Wcdt=0.Substituting equations (6), (10), and (14) into equation (15) and ignoring the nonlinear term in the equation, the vibration differential equation of the uniformly rotating (θ˙=ω0,θ¨=0) fractional derivative viscoelastic beam can be obtained:(16)−EIw,xxxx+ρA−w¨+ω02w−23ρIω02w,xx+23ρIw¨,xx+ρAω02aL−xw,xx−ηIDαw,xxxx+12ρAω02L2−x2w,xx−ρAω02a+xw,x=0.The boundary conditions of the fractional derivative viscoelastic beam are as follows:x=0:(17a)w=0,(17b)w,x=0.x=L:(17c)EIw,xx+ηIDαw,xx=0,(17d)EI∂∂x∂2w∂x2−23ρI∂3w∂x∂t2+23ρIω02∂w∂x+ηI∂∂xDαw,xx=0.


### 3.3. Dimensional Normalization of the Equation

To facilitate subsequent calculations, the following dimensionless quantities have been introduced:(18)ξ=xL,W=wL,r0=aL,r=RL,t˜=tL2ρA/EI,ω1=ωL02ρAEI,g=ηEL2ρA/EI,where r0 is the hub radius-beam length ratio, r is slenderness ratio, ω1 is dimensionless hub speed, and g is viscoelastic ratio.Substituting the above dimensionless quantity into equation (16)(19)∂4W∂ξ4+∂2W∂t˜2−ω12r01−ξ∂2W∂ξ2+16ω12r2∂2W∂ξ2−12ω121−ξ2∂2W∂ξ2−ω12W+g∂α∂t˜α∂4W∂ξ4−16r2∂4W∂ξ2∂t˜2+ω12r0+ξ∂W∂ξ=0.The solution of equation (19) can be set as Wξ,t˜=Wξexpωft˜, and the differential equation of mode shape can be written as follows:(20)W,ξξξξξ+ωf2Wξ−ω12Wξ+16ω12r2W,ξξξ−ω12r01−ξ+121−ξ2W,ξξξ+gωfαW,ξξξξξ+ω12r0+ξW,ξξ−16ωf2r2W,ξξξ=0,where ωf is a dimensionless complex frequency.Substituting the dimensionless quantity into equations (17a) and –(17d), respectively, the dimensionless boundary conditions of the fractional derivative viscoelastic beam can be expressed as follows:(21a)W0=0,(21b)W,ξ0=0,(21c)W,ξξ1+gωfαW,ξξ1=0,(21d)W,ξξξ1−16r2ωf2W,ξ1+16r2ω12W,ξ1+gωfαW,ξξξ1=0.


### 3.4. Differential Quadrature Method

When the differential quadrature method is used to solve the equation, most of the calculation accuracy will be determined by the selection of the weight function. Therefore, when using this method for calculation, the node distribution must be selected reasonably. For the equations of the specific problem to be solved, the properties of the solution region are different, often using nonuniform node partitioning, and more mesh points can be selected for regions close to the boundary where the physical properties change more prominently. And for the boundary region, the small variableδ is used for the adjacency processing.The values of the nodes by nonuniform sampling method are as follows:(22)ξ1=0,ξ2=δ,ξi=121−cosi−2N−3π,i=3,4,…,N−2,ξN−1=1−δ,ξN=1,where set δ=1×10−6 and N=13 in subsequent calculation.Then the discretization of equation (20) can be written as follows:(23)∑j=1NAij4Wj+16ω12r2∑j=1NAij2Wj+gωfα∑j=1NAij4Wj−ω12r01−ξi+121−ξi2∑j=1NAij2Wj−ω12Wi+ωf2Wi+ω12r0+ξi∑j=1NAij1Wj−16ωf2r2∑j=1NAij2Wj=0,where Aij1, Aij2, Aij3, and Aij4 are first-order, second-order, third-order, and fourth-order weight coefficients, respectively,(24)Aij1=∏k=1k≠i,jNxi−xk∏k=1k≠jNxj−xk,i,j=1,2,…,N;i≠j,∑k=1k≠iN1xi−xk,i,j=1,2,…,N;i=j,Aij2=∑k=1NAik1Akj1,Aij3=∑k=1NAik1Akj2,i,j=1,2,…,N,Aij4=∑k=1NAik2Akj2.The discretization of the boundary conditions is(25a)W1=0,(25b)∑j=1NA2j1Wj=0,(25c)1+gωfα∑j=1NAN−1,j2Wj=0,(25d)∑j=1NANj3Wj−16r2ωf2∑j=1NANj1Wj+16r2ω12∑j=1NANj1Wj+gωfα∑j=1NANj3Wj=0.Equations (23), (25a) and –(25d) can be written in matrix form:(26)ωf2M+ωfαC+KW=0,where M, C, and K are all N−1 square matrix and W=W2W3⋯WN−1WNT is column matrix.Thus, the generalized characteristic equation corresponding to equation (26) can be obtained as follows:(27)ωf2M+ωfαC+K=0,where M, C and K contain the hub radius-beam length ratio, the slenderness ratio, the viscoelastic ratio, the dimensionless hub speed, fractional order, and dimensionless complex frequency of the system.


## 3.1. Deformation Field Description

Figure1(a) shows that the schematic diagram of the fractional derivative viscoelastic rotating hollow beam, which has length L of the beam, radius a, and rotating angular velocity ω of the hub. The oxy and OXY are the floating coordinate system and the inertial coordinate system, respectively. The θ represents rotation angle of the viscoelastic beam in a wide range of motion, that is, the rotation angle of the oxy is relative to the OXY. The ox-axis coincides with the axis of the viscoelastic hollow circular section beam before deformed.Figure 1
Dynamic model of fractional derivative viscoelastic rotating hollow beam. (a) Fractional derivative viscoelastic rotating hollow beam. (b) Viscoelastic rotating hollow beam profile.
(a)
(b)In Figure1(b), the cross-section of the beam is a hollow section with wall thickness h and Rn is the radius of the wall thickness for any cross-section at the middle line. In the n−s coordinate system, n is the normal direction and s is the tangential direction.In Figure2, ra is a radius vector of the hub, xp is the vector diameter of any point P on the axis of a fractional derivative viscoelastic beam relative to the point o in the floating coordinate system when it is not deformed, and up is a displacement vector of the point P.Figure 2
Displacement field at any point on the beam axis.The position vector of pointP relative to the point O can be written as follows:(3)rp=ra+Aθxp+up,where(4)ra=acosθasinθT,xp=x0T,Aθ=cosθ−sinθsinθcosθ,up=upxupy=u+uc+ubw=u−12∫0xw,x2dx−yw,xw,where Aθ, u, and w are the transformation matrices of the oxy relative to the OXY, the axial displacement, and the lateral bending deformation displacement on the beam axis, respectively; ub=−yw,x is the axial displacement caused by the rotation of the beam cross-section, uc=−1/2∫0xw,x2dx is the secondary deformation coupling of the beam, that is, the axial expansion amount caused by the lateral bending deformation of the fractional derivative viscoelastic rotating hollow beam, and “,x” represents the first derivative of x.Taking the first derivative of time for equation (3), we get the velocity vector of point P as follows:(5)r˙p=r˙a+θ˙×Αθxp+up+Αθu˙p=−a+x+uc+ubθ˙sinθ−wθ˙cosθ+u˙c+u˙bcosθ−w˙sinθa+x+uc+ubθ˙cosθ−wθ˙sinθ+u˙c+u˙bsinθ+w˙cosθ,where θ˙=0−θ˙θ˙0 is the antisymmetric matrix of angular velocity θ˙ and “⋅” represents the first derivative of t.


## 3.2. Differential Equations of Motion

The total kinetic energy of the system can be seen as the sum of the kinetic energy of the hub and the kinetic energy of the fractional derivative viscoelastic rotating beam, namely(6)Ts=TH+T=12JHθ˙2+12ρ∫Vr˙pTr˙pdV=12JHθ˙2+13ρIθ˙2∫0Lw,x2dx+13ρI∫0Lw˙,x2dx+12∫0LρAw˙2+u˙c2+θ˙2w2+a+x+uc2dx+∫0LρAθ˙−u˙cw+w˙a+x+ucdx,where JH is the inertia moment of the hub with respect to the central axis.The strain energy of the system can be divided into two parts, one part is potential energy and the other part is dissipative energy.The potential energy can be expressed as follows:(7)V0=∫V12σxεxdV=∫V12Eεx2dV,where σx, εx, and V are the normal stress, normal strains of the beam, and the volume of the viscoelastic beam, respectively.In equation (4), not considering the axial displacement and ignoring nonlinear term, the normal strain can be obtained by the relationship between strain and displacement, namely,(8)εx=upx,x=−12∫0xw,x2dx−yw,x,x=−ys−nz,sw,xx−12w,x2.According to Figure1(b), the geometric relations can be obtained:(9)dzds=−sinϕ=−sinsRn,dyds=cosϕ=cossRn.Then, the potential energy can be expressed as follows:(10)Vo=∫0L∫A12Eεx2tdAdx=∫0L∫−h/2h/2∫0πRn12ERn+nsinsRnw,xx+12w,x22Rn+nRndsdndx=12EI∫0Lw,xx2dx+18EA∫0Lw,x4dx,where I and A are the inertia moment and the cross-sectional area of the beam, respectively,(11)I=164π2Rn+h4−2Rn−h4=14π4Rn3h+Rnh3,A=πRn+h22−Rn−h22=2πRnh.The expression of the system’s dissipative energy is(12)Wc=∫V∫0εxtηDαεxtdεxtdV.From equation (2), the following equation can be obtained(13)Dαεxt=εx0Γ1−αt−α+1Γ1−α∫0tε˙xτt−ταdτ.Substituting equation (12) into equation (13), the following equation can be obtained:(14)Wc=∫VηDαεxtεxtdV=∫V∫0εxtηεx0Γ1−αt−αdεxtdV+∫V∫0εxtη1Γ1−α∫0tε˙xτt−ταdτdεxxdV=ηI∫0Lw,xxDαw,xxdx+14ηA∫0Lw,x2Dαw,x2dx.In this paper, based on Hamiltonian variational principle, the differential equation of fractional derivative viscoelastic rotating beam is derived. Hamilton’s variational principle takes the form of(15)δ∫t1t2Ldt=δ∫t1t2Ts−V0−Wcdt=0.Substituting equations (6), (10), and (14) into equation (15) and ignoring the nonlinear term in the equation, the vibration differential equation of the uniformly rotating (θ˙=ω0,θ¨=0) fractional derivative viscoelastic beam can be obtained:(16)−EIw,xxxx+ρA−w¨+ω02w−23ρIω02w,xx+23ρIw¨,xx+ρAω02aL−xw,xx−ηIDαw,xxxx+12ρAω02L2−x2w,xx−ρAω02a+xw,x=0.The boundary conditions of the fractional derivative viscoelastic beam are as follows:x=0:(17a)w=0,(17b)w,x=0.x=L:(17c)EIw,xx+ηIDαw,xx=0,(17d)EI∂∂x∂2w∂x2−23ρI∂3w∂x∂t2+23ρIω02∂w∂x+ηI∂∂xDαw,xx=0.


## 3.3. Dimensional Normalization of the Equation

To facilitate subsequent calculations, the following dimensionless quantities have been introduced:(18)ξ=xL,W=wL,r0=aL,r=RL,t˜=tL2ρA/EI,ω1=ωL02ρAEI,g=ηEL2ρA/EI,where r0 is the hub radius-beam length ratio, r is slenderness ratio, ω1 is dimensionless hub speed, and g is viscoelastic ratio.Substituting the above dimensionless quantity into equation (16)(19)∂4W∂ξ4+∂2W∂t˜2−ω12r01−ξ∂2W∂ξ2+16ω12r2∂2W∂ξ2−12ω121−ξ2∂2W∂ξ2−ω12W+g∂α∂t˜α∂4W∂ξ4−16r2∂4W∂ξ2∂t˜2+ω12r0+ξ∂W∂ξ=0.The solution of equation (19) can be set as Wξ,t˜=Wξexpωft˜, and the differential equation of mode shape can be written as follows:(20)W,ξξξξξ+ωf2Wξ−ω12Wξ+16ω12r2W,ξξξ−ω12r01−ξ+121−ξ2W,ξξξ+gωfαW,ξξξξξ+ω12r0+ξW,ξξ−16ωf2r2W,ξξξ=0,where ωf is a dimensionless complex frequency.Substituting the dimensionless quantity into equations (17a) and –(17d), respectively, the dimensionless boundary conditions of the fractional derivative viscoelastic beam can be expressed as follows:(21a)W0=0,(21b)W,ξ0=0,(21c)W,ξξ1+gωfαW,ξξ1=0,(21d)W,ξξξ1−16r2ωf2W,ξ1+16r2ω12W,ξ1+gωfαW,ξξξ1=0.


## 3.4. Differential Quadrature Method

When the differential quadrature method is used to solve the equation, most of the calculation accuracy will be determined by the selection of the weight function. Therefore, when using this method for calculation, the node distribution must be selected reasonably. For the equations of the specific problem to be solved, the properties of the solution region are different, often using nonuniform node partitioning, and more mesh points can be selected for regions close to the boundary where the physical properties change more prominently. And for the boundary region, the small variableδ is used for the adjacency processing.The values of the nodes by nonuniform sampling method are as follows:(22)ξ1=0,ξ2=δ,ξi=121−cosi−2N−3π,i=3,4,…,N−2,ξN−1=1−δ,ξN=1,where set δ=1×10−6 and N=13 in subsequent calculation.Then the discretization of equation (20) can be written as follows:(23)∑j=1NAij4Wj+16ω12r2∑j=1NAij2Wj+gωfα∑j=1NAij4Wj−ω12r01−ξi+121−ξi2∑j=1NAij2Wj−ω12Wi+ωf2Wi+ω12r0+ξi∑j=1NAij1Wj−16ωf2r2∑j=1NAij2Wj=0,where Aij1, Aij2, Aij3, and Aij4 are first-order, second-order, third-order, and fourth-order weight coefficients, respectively,(24)Aij1=∏k=1k≠i,jNxi−xk∏k=1k≠jNxj−xk,i,j=1,2,…,N;i≠j,∑k=1k≠iN1xi−xk,i,j=1,2,…,N;i=j,Aij2=∑k=1NAik1Akj1,Aij3=∑k=1NAik1Akj2,i,j=1,2,…,N,Aij4=∑k=1NAik2Akj2.The discretization of the boundary conditions is(25a)W1=0,(25b)∑j=1NA2j1Wj=0,(25c)1+gωfα∑j=1NAN−1,j2Wj=0,(25d)∑j=1NANj3Wj−16r2ωf2∑j=1NANj1Wj+16r2ω12∑j=1NANj1Wj+gωfα∑j=1NANj3Wj=0.Equations (23), (25a) and –(25d) can be written in matrix form:(26)ωf2M+ωfαC+KW=0,where M, C, and K are all N−1 square matrix and W=W2W3⋯WN−1WNT is column matrix.Thus, the generalized characteristic equation corresponding to equation (26) can be obtained as follows:(27)ωf2M+ωfαC+K=0,where M, C and K contain the hub radius-beam length ratio, the slenderness ratio, the viscoelastic ratio, the dimensionless hub speed, fractional order, and dimensionless complex frequency of the system.


## 4. Numerical Results and Analyses

### 4.1. Comparison and Verification

The correctness and feasibility of the differential quadrature method used in this paper are analyzed by comparing the numerical results with reference [7]. In the case when g=0, α=0, r=0, and r0=0.125, set the dimensionless hub speed at 0, 5, and 10, the result of the compared first three-order natural frequency with reference [7] is shown in Table 1:Table 1
Comparison of the first three dimensionless natural frequencies of elastic cantilever rotating beams atg=0, α=0, r=0, and r0=0.125.
ω
1
0
5
10
Present
Reference [7]
Present
Reference [7]
Present
Reference [7]
1st
3.515
3.507
6.302
6.425
10.953
11.159
2nd
22.029
21.647
24.951
24.973
32.122
32.965
3rd
60.673
59.207
65.011
63.539
73.672
71.497It can be seen from the Table1 that when the dimensionless hub rotational speed ω1 is 0, 5, and 10 respectively, the dimensionless natural frequency of the system is basically the same as that of reference [7], and all errors are below 2%.


### 4.2. Numerical Calculation and Analysis

Considering the influence of moment of inertia, the effects of slenderness ratior, the hub radius-beam length ratio r0, viscoelastic ratio g, dimensionless hub speed ω1, and fractional order α on the first three-order dimensionless Imωf (imaginary part of complex frequency) and Reωf (real part of complex frequency) of the fractional derivative rotating viscoelastic beam are analyzed.Figure3 plots the curves of complex frequency of the system with the hub radius-beam length ratio at r=1/30, g=0.025, ω1=10, and α=1/2. And we can see that the complex frequency of the system vary linearly with r0, and r0 is positively correlated with the Imωf, but negatively correlated with the Reωf. With the increase of the order, the effect of r0 on the Reωf is more obvious.Figure 3
Variation curves of complex frequency of the system with the hub radius-beam length ratio atr=1/30, g=0.025, ω1=10, and α=1/2.
(a)
(b)Figure4 plots the curves of complex frequency of the system with slenderness ratio at r0=0.25, g=0.025, ω1=10, and α=1/2.Figure 4
Variation curves of complex frequency of the system with slenderness ratio atr0=0.25, g=0.025, ω1=10, and α=1/2.
(a)
(b)And we can see that when the slenderness ratio is between 0 and 0.08, slenderness ratio has less influence on the complex frequency of the system, and when slenderness ratio is between 0.08 and 0.10, theImωf of the system is negatively correlated with slenderness ratio, and the Reωf increases first and then decreases. With the increase of the order of the system, the peak value of the Reωf moves forward.Figure5 plots the curves of complex frequency of the system with viscoelastic ratio at r0=0.25, r=1/30, ω1=10, and α=1/2.Figure 5
Variation curves of complex frequency of the system with viscoelastic ratio atr0=0.25, r=1/30, ω1=10, and α=1/2.
(a)
(b)And we can see that the viscoelastic ratio has slight effect on theImωf of the system, and with the increase of the viscoelastic ratio, the Imωf hardly changes. Meanwhile, the viscoelastic ratio has a great influence on the Reωf, and there is a positive correlation between them, and with the increase of order, the Reωf increases more obviously. Therefore, the viscoelastic ratio can be used to describe the damping characteristics of materials.Figure6 plots the curve of complex frequency of the system with dimensionless hub speed at different fractional order α=1/10, α=1/5, and α=1/2 at r0=0.25, r=1/30, and g=0.025. And we can see that the dimensionless hub speed has a significant influence on the complex frequency of the system, and as the speed of the dimensionless hub increases, the Imωf increases, the Reωf decreases, and the attenuation of higher-order Reωf becomes more obvious.Figure 6
Variation curves of complex frequency of the systems at dimensionless hub speed with different fractional order atr0=0.25, r=1/30, and g=0.025.
(a)
(b)Meanwhile, fractional order has little effect on theImωf of the system, but the influence on the Reωf is great, and the influence of fractional order on the high-order Reωf is obvious. When α=1/2, the attenuation rate of the Reωf is faster than α=1/5. Thus, viscoelastic materials with larger fractional order will be more conducive to the stability of the system.


### 4.3. Response of Viscoelastic Rotating Beams with Kelvin–Voigt Fractional Derivatives

For the system, set the number of mesh pointsN=13, small variables δ=1×10−6, and nonuniform meshing of Kelvin–Voigt viscoelastic material beam with fractional derivative, when ξ3=0.02447, ξ7=0.5, ξ11=0.97553, and ξ13=1, and the response curve of different points on the viscoelastic beam is as follows Figure 7:Figure 7
Vibration response curve of fractional derivative Kelvin–Voigt viscoelastic beam atr=1/30, r0=0.25, g=0.025, ω1=10, and α=1/2. (a) ξ3=0.02447, (b) ξ7=0.5, (c) ξ11=0.97553, and (d) ξ13=1.
(a)
(b)
(c)
(d)Figure7 plots the curves of response of the Kelvin–Voigt viscoelastic beam with fractional derivative when r0=0.25, r=1/20, g=0.025, ω1=10, and α=1/2. And we can see that the response at different points of the Kelvin–Voigt viscoelastic material beam.Asξi increases, the corresponding response amplitude also increases, which indicates that the closer to the free end, the greater the response amplitude. Moreover, the amplitude of the system’s response decreases with time, which is due to the damping effect in the material.


## 4.1. Comparison and Verification

The correctness and feasibility of the differential quadrature method used in this paper are analyzed by comparing the numerical results with reference [7]. In the case when g=0, α=0, r=0, and r0=0.125, set the dimensionless hub speed at 0, 5, and 10, the result of the compared first three-order natural frequency with reference [7] is shown in Table 1:Table 1
Comparison of the first three dimensionless natural frequencies of elastic cantilever rotating beams atg=0, α=0, r=0, and r0=0.125.
ω
1
0
5
10
Present
Reference [7]
Present
Reference [7]
Present
Reference [7]
1st
3.515
3.507
6.302
6.425
10.953
11.159
2nd
22.029
21.647
24.951
24.973
32.122
32.965
3rd
60.673
59.207
65.011
63.539
73.672
71.497It can be seen from the Table1 that when the dimensionless hub rotational speed ω1 is 0, 5, and 10 respectively, the dimensionless natural frequency of the system is basically the same as that of reference [7], and all errors are below 2%.


## 4.2. Numerical Calculation and Analysis

Considering the influence of moment of inertia, the effects of slenderness ratior, the hub radius-beam length ratio r0, viscoelastic ratio g, dimensionless hub speed ω1, and fractional order α on the first three-order dimensionless Imωf (imaginary part of complex frequency) and Reωf (real part of complex frequency) of the fractional derivative rotating viscoelastic beam are analyzed.Figure3 plots the curves of complex frequency of the system with the hub radius-beam length ratio at r=1/30, g=0.025, ω1=10, and α=1/2. And we can see that the complex frequency of the system vary linearly with r0, and r0 is positively correlated with the Imωf, but negatively correlated with the Reωf. With the increase of the order, the effect of r0 on the Reωf is more obvious.Figure 3
Variation curves of complex frequency of the system with the hub radius-beam length ratio atr=1/30, g=0.025, ω1=10, and α=1/2.
(a)
(b)Figure4 plots the curves of complex frequency of the system with slenderness ratio at r0=0.25, g=0.025, ω1=10, and α=1/2.Figure 4
Variation curves of complex frequency of the system with slenderness ratio atr0=0.25, g=0.025, ω1=10, and α=1/2.
(a)
(b)And we can see that when the slenderness ratio is between 0 and 0.08, slenderness ratio has less influence on the complex frequency of the system, and when slenderness ratio is between 0.08 and 0.10, theImωf of the system is negatively correlated with slenderness ratio, and the Reωf increases first and then decreases. With the increase of the order of the system, the peak value of the Reωf moves forward.Figure5 plots the curves of complex frequency of the system with viscoelastic ratio at r0=0.25, r=1/30, ω1=10, and α=1/2.Figure 5
Variation curves of complex frequency of the system with viscoelastic ratio atr0=0.25, r=1/30, ω1=10, and α=1/2.
(a)
(b)And we can see that the viscoelastic ratio has slight effect on theImωf of the system, and with the increase of the viscoelastic ratio, the Imωf hardly changes. Meanwhile, the viscoelastic ratio has a great influence on the Reωf, and there is a positive correlation between them, and with the increase of order, the Reωf increases more obviously. Therefore, the viscoelastic ratio can be used to describe the damping characteristics of materials.Figure6 plots the curve of complex frequency of the system with dimensionless hub speed at different fractional order α=1/10, α=1/5, and α=1/2 at r0=0.25, r=1/30, and g=0.025. And we can see that the dimensionless hub speed has a significant influence on the complex frequency of the system, and as the speed of the dimensionless hub increases, the Imωf increases, the Reωf decreases, and the attenuation of higher-order Reωf becomes more obvious.Figure 6
Variation curves of complex frequency of the systems at dimensionless hub speed with different fractional order atr0=0.25, r=1/30, and g=0.025.
(a)
(b)Meanwhile, fractional order has little effect on theImωf of the system, but the influence on the Reωf is great, and the influence of fractional order on the high-order Reωf is obvious. When α=1/2, the attenuation rate of the Reωf is faster than α=1/5. Thus, viscoelastic materials with larger fractional order will be more conducive to the stability of the system.


## 4.3. Response of Viscoelastic Rotating Beams with Kelvin–Voigt Fractional Derivatives

For the system, set the number of mesh pointsN=13, small variables δ=1×10−6, and nonuniform meshing of Kelvin–Voigt viscoelastic material beam with fractional derivative, when ξ3=0.02447, ξ7=0.5, ξ11=0.97553, and ξ13=1, and the response curve of different points on the viscoelastic beam is as follows Figure 7:Figure 7
Vibration response curve of fractional derivative Kelvin–Voigt viscoelastic beam atr=1/30, r0=0.25, g=0.025, ω1=10, and α=1/2. (a) ξ3=0.02447, (b) ξ7=0.5, (c) ξ11=0.97553, and (d) ξ13=1.
(a)
(b)
(c)
(d)Figure7 plots the curves of response of the Kelvin–Voigt viscoelastic beam with fractional derivative when r0=0.25, r=1/20, g=0.025, ω1=10, and α=1/2. And we can see that the response at different points of the Kelvin–Voigt viscoelastic material beam.Asξi increases, the corresponding response amplitude also increases, which indicates that the closer to the free end, the greater the response amplitude. Moreover, the amplitude of the system’s response decreases with time, which is due to the damping effect in the material.


## 5. Conclusion

Based on Hamilton principle and Euler–Bernoulli beam theory, the influence of slenderness ratio, the hub radius-beam length ratio, dimensionless hub speed, viscoelastic ratio, and fractional order on the vibration characteristics of viscoelastic rotating beams with fractional derivatives is investigated by the differential quadrature method. The calculation results in this paper have little error with the existing literature, which shows the feasibility and correctness of the differential quadrature method. Compared with other methods, the solving process of the differential quadrature method is more intuitive and easier to program.When the dimensionless hub speed is constant, the dimensionlessImωf and Reωf of the system both increases slightly with the hub radius-beam length ratio increase. For a fractional derivative viscoelastic rotating hollow beam at a constant velocity, when slenderness ratio is about 0.08, the Imωf and Reωf have a sudden change; the former has been attenuating, while the latter has increased first and then decreased rapidly. The viscoelastic ratio has little effect on the Imωf of the system, but it has a great influence on the Reωf, and at high order, the Reωf increases rapidly.With the increase of dimensionless hub speed, theImωf of the system decreases, while the Reωf increases, and the attenuation of higher-order Reωf is more significant. Fractional order has little effect on the Imωf of the system, but it has great influence on the Reωf, and the effect on the higher-order Reωf is more obvious.Through the study of response at different points on a fractional derivative Kelvin–Voigt viscoelastic beam, we can see that the response at different points on a viscoelastic beam varies with time in a wave function, and the closer to the fixed end, the smaller the amplitude of response. When near the free end, the amplitude of response is larger but the amplitude of the system's response decreases with time.



---

*Source: 5715694-2019-10-16.xml*