Jtam-A4.dvi JOURNAL OF THEORETICAL AND APPLIED MECHANICS 54, 4, pp. 1433-1445, Warsaw 2016 DOI: 10.15632/jtam-pl.54.4.1433 DYNAMIC RESPONSE OF A SIMPLY SUPPORTED VISCOELASTIC BEAM OF A FRACTIONAL DERIVATIVE TYPE TO A MOVING FORCE LOAD Jan Freundlich Warsaw University of Technology, Faculty of Automotive and Construction Machinery Engineering, Warszawa, Poland e-mail: jfr@simr.pw.edu.pl In the paper, the dynamic response of a simply supported viscoelastic beam of the fractio- nal derivative type to a moving force load is studied. The Bernoulli-Euler beam with the fractional derivative viscoelastic Kelvin-Voigt material model is considered. The Riemann- Liouville fractional derivative of the order 0<α¬ 1 is used. The forced-vibration solution of the beam is determined using the mode superposition method. A convolution integral of fractional Green’s function and forcing function is used to achieve the beam response. Green’s function is formulated by two terms. The first term describes damped vibrations around the drifting equilibrium position, while the second term describes the drift of the equilibrium position. The solution is obtained analytically whereas dynamic responses are calculated numerically. A comparison between the results obtained using the fractional and integer viscoelastic material models is performed. Next, the effects of the order of the frac- tional derivative and velocity of the moving force on the dynamic response of the beam are studied. In the analysed system, the effect of the termdescribing the drift of the equilibrium position on the beamdeflection is negligible in comparisonwith the first term and therefore can be omitted. The calculated responses of the beam with the fractional material model are similar to those presented in works of other authors. Keywords: fractional viscoelasticity, moving loads, beam vibrations, transient dynamic analysis 1. Introduction The problem of predicting the dynamic response of a structure resulting from the passage of moving loads often occurs in engineering analysis. This problem has many applications in en- gineering analysis, some of these applications are vibrations that occur in bridges and railroad tracks, cranes, machine elements, weapon barrels, rails, rail and vehicle body parts. The studies show that the transversal deflection and stresses due to moving loads are considerably higher than those observed with stationary loads. The dynamic behaviour of a beam subjected tomo- ving loads or moving masses have been studied extensively in connection with the design of railway tracks and bridges as well as machining processes. A broad overview of the analytical and numerical methods deals with these issues is provided in a book by Fryba (1972). A newer survey of these methods can be found e.g. in a book by Bajer andDyniewicz (2012). The book also presents a comprehensive overview of the numerical methods including the finite element method (FEM) and the space-time element method (STEM). Additionally, both these books contain a broad bibliography devoted to the dynamic analysis of structures subjected tomoving loads. Furthermore, damping effects on structure dynamics are considered in many works devoted to the mentioned above subject (e.g Fryba, 1972). Therefore, selection of an appropriate visco- elasticmaterialmodel is an essential element in studying the dynamical behaviour ofmechanical structures. The selected model should as accurately as possible describe viscoelastic damping 1434 J. Freundlich over a wide interval of frequencies. A large number of engineering materials show a weak fre- quency dependence of their damping properties within a broad frequency range (Bagley and Torvik, 1983a,b; Caputo andMainardi, 1971a; Clough andPenzien, 1993). This weak frequency dependence is difficult to describe using classical viscoelasticmodels, based on integer order rate laws. In the recent years, the fractional derivatives have been frequently utilized to describe dynamic characteristics of viscoelastic materials and damping elements (Podlubny, 1999; Ma- inardi, 2009; Rossikhin and Shitikova, 2010; Hedrih (Stevanović) 2014; Hedrih (Stevanović) and Machado, 2015). By introducing fractional-order derivatives instead of integer-order derivatives in the constitutive relations, the number of parameters required to accurately describe the dyna- mic properties can be significantly reduced (see, e.g., Caputo andMainardi, 1971b; Bagley and Torvik, 1983b; Enelund and Olsson, 1999). A review of publications dedicated to application of fractional calculus in dynamics of solids can be found in the paper by Rossikhin and Shiti- kova (2010). Moreover, a fractional damping is utilized in dynamical analysis of beams under moving loads (e.g. Abu-Mallouh et al., 2012; Alkhaldi et al., 2013). In these papers, taken into consideration are only loadsmoving at a constant speed. The solutions presented in these works are obtained with the help of Mittag-Lefler functions and infinite series. In numerous studies in which fractional viscoelastic models are employed, the dynamic transient analysis is limited to the response induced by initial conditions. However, structures subjected tomoving loadswith a time varying velocity are often encountered in engineering practice. The number of publications considering the dynamic behaviour of a beamwith a fractional viscoelstic model under moving loads is rather limited. Numerous research works related to structure vibrations described by fractional models present transient solutions expressed as series expansions (Atanackovic and Stankovic, 2002; Bagley and Torvik, 1986; Caputo, 1974; Hedrih (Stevanović) and Filipovski, 2002; Podlubny, 1999) which are not very useful for practical applications because of their slow convergence rate (Podlubny, 1999; Rossikhin and Shitikova, 2010). Thus, the aim of this work is demonstration of an approach utilizing fractional Green’s function in dynamical analysis of a simply-supported beam subjected to a moving force with a constant velocity and constant acceleration. Green’s function is evaluated using a closed contour of integration in conjunction with the residue theorem (Caputo and Mainardi, 1971b; Bagley and Torvik, 1983b, 1986; Ros- sikhin and Shitikova, 1997; Beyer andKempfle, 1995). The proposed approach is the expansion of methods employing fractional calculus in the modelling of damping properties of dynamic systems. 2. Problem formulation The equation of motion of the examined beam is derived on the assumption of the Bernoulli- -Euler theory, neglecting rotary inertia and shear deformation. Furthermore, it is assumed that the beam is homogeneous with constant cross-section along its length, and the neutral axis of the beambending is one of the principal axes of inertia of the normal cross section of the beam. Beamoscillations are assumed in thexy plane (Fig. 1).Moreover, it is assumed that the internal dissipation of mechanical energy in the beammaterial is described by a differential equation of the fractional order. Therefore, the stress-strain constitutive relation of the beammaterial is in the following form (e.g. Bagley and Torvik, 1983a) σ=Eε(t)+E′γD γ t [ε(t)] =E ( ε+µγ dγε(t) dtγ ) (2.1) where: µγ = E ′ γ/E, E – Young’s modulus of the beam material, E ′ γ – damping coefficient, D γ t [·] – fractional order differential operator of the γ-th derivative with respect to time t in the following form (Miller and Ross, 1993; Podlubny, 1999) Dynamic response of a simply supported viscoelastic beam... 1435 D γ t f(t)≡ dγ dtγ f(t)≡ 1 Γ(1−γ) d dt t∫ 0 f(τ) dτ (t− τ)γ (2.2) whereΓ(1−γ) is the Euler gamma function and t­ 0 (Miller andRoss, 1993; Podlubny, 1999). For many real materials, the fractional derivative order is commonly assumed to be in the interval 0 < γ < 1 (Caputo and Mainardi, 1971a; Bagley and Torvik, 1983b; Enelund and Olsson, 1999). The constitutive relation with γ = 1 represents Kelvin-Voight material with internal linear dissipation ofmechanical energy (e.g.Mainardi, 2009; Hedrih (Stevanović), 2014) Fig. 1. Scheme of the analysed system The beam is assumed to be subjected to a forcemovingwith a constant velocity or accelera- tion from the left to the right edge of the beam, according with the sense of the x axis (Fig. 1). Based on the above assumptions, the governing equation of the analysed beam is obtained as below EJ (∂4y(x,t) ∂x4 +µγ dγ dtγ (∂4y(x,t) ∂x4 )) +Aρ ∂2y(x,t) ∂t2 =Fδ(x−g(t)) (2.3) where A is the cross-section area of the beam, J – axial moment of inertia of the beam cross- -section respect to the neutral axis of the beambending, ρ –material mass density of the beam, y(x,t) – transversal displacement of the neutral beam axis (Fig. 1), t – time, x – longitudinal coordinate,F – intensity of the externalmoving and concentrated force, δ –Dirac delta function, g(t) – function determining the force location, g(t) = vt or g(t) = εt2/2 for a forcemovingwith a constant velocity or constant acceleration, respectively. Furthermore, this functionmust satisfy the following condition 0¬ g(t)¬ l. For the simply-supportedbeam, the deflection and bendingmoment at both beam ends have to be zero, thus the boundary conditions are as below y(0, t) = 0 y(l,t)= 0 Mb(0, t) = 0 Mb(l,t)= 0 (2.4) where l is the length of the beam,Mb – bendingmoment. The bendingmoment is calculated as below Mb(x,t)= ∫ A σy dA= ∫ A {Eε(t)+E′γD γ t [ε(t)]}y dA=−EJ {d2y(x,t) dx2 +µγD γ t [d2y(x,t) dx2 ]} (2.5) therefore d2y(0, t) dx2 = d2y(l,t) dx2 =0 D γ t [d2y(0, t) dx2 ] =D γ t [d2y(l,t) dx2 ] =0 (2.6) Moreover, it is assumed that the beam is initially at rest and the initial bending moments are equal to 0. 1436 J. Freundlich The solution to Eq. (2.3) is obtained utilizing themode superposition principle (e.g. Kaliski, 1966; Clough and Penzien, 1993; Rao, 2004). The eigenfunctions for a simply supported beam are given by Yn(x)= sin nπx l n=1,2,3, . . . (2.7) Then, the forced-vibration solution of a beam can be expressed as y(x,t)= ∞∑ n=1 ξn(t)Yn(x)= ∞∑ n=1 ξn(t)sin nπx l (2.8) and the corresponding derivatives are evaluated below ∂4y(x,t) ∂x4 = ∞∑ n=1 Y IVn (x)ξn(t) dγ dt (∂4y(x,t) ∂x4 ) = ∞∑ n=1 Y IVn (x)D γ t [ξ(t)] (2.9) where D γ t [ξn(t)] = dγξ(t) dtγ = 1 Γ(1−γ) d dt t∫ 0 ξn(t) dτ (t− τ)γ and Y IVn (x)= d4Y (x) dx4 ∂2y(x,t) ∂t2 = ∞∑ n=1 Yn(x)ξ̈n(t) (2.10) where ξ̈(t)= d2ξ(t) dt2 (2.11) Substituting Eqs. (2.9) and (2.10) into (2.3), we get ∞∑ n=1 ( Y IVn (x)ξn(t)+µγY IV n (x)D γ t [ξn(t)]+a 2 bY (x)ξ̈n(t) ) = F EJ δ(x−g(t)) (2.12) Next, integrating over x from 0 to 1, and using the orthogonality property of eigenfunctions, after somemathematical transformations, the following relationships are obtained ξ̈n(t)+µγω 2 nD γ t [ξn(t)]+ω 2 nξn(t)= f0 sin nπg(t) l (2.13) where f0 = 2F m ωn = (πn l )2 √ EJ ρA m – mass of the beam. Assuming zero initial conditions ξ(t) ∣∣∣ t=0 =0 d[ξ(t)] dt ∣∣∣∣∣ t=0 =0 dγ−1[ξ(t)] dtγ−1 ∣∣∣∣∣ t=0 =0 the solution of equations (2.13) can be expressed as ξn(t)= f0 t∫ 0 Gn(t− τ)sin nπg(τ) l dτ (2.14) Dynamic response of a simply supported viscoelastic beam... 1437 whereGn(t) is Green’s function corresponding to Eq. (2.13) (Miller and Ross, 1993; Podlubmy, 1999; Rossikhin and Shitikova, 1997, 2010). ThisGreen’s function consists of two terms, namely Gn(t)=K1n(t)+K2n(t) (2.15) The first term K1n (Eq. (2.15)) represents damped vibrations around the drifting equilibrium position, while the second term K2n describes the drift of the equilibrium position (Beyer and Kempfle, 1995; Rossikhin and Shitikova, 1997). The term K1n could be calculated from the formula given by Beyer and Kempfle (1995) K1n(t)=αne −σnt sin(Ωnt+φn) (2.16) where αn = 2 √ µ2 k +ν2 φn =arctan µk ν µk =Re(W ′ (p,1,2)) ν =Im(W ′ (p,1,2)) and W(p) = p2 + µγω 2 np γ + ω2n is the characteristic polynomial of equations (2.13), W ′(p) = 2p + γµγω 2 np γ−1 – derivative of the characteristic polynomial with respect p, p1,2 =−σn± iΩn – conjugate complex roots of the characteristic polynomialW(p). The termK2n could be calculated using the formula (Beyer andKempfle, 1995) K2n(t)= µγω 2 n sin(πγ) π ∞∫ 0 rγe−rt dr [r2+µγω 2 nr γ cos(πγ)+ω2n] 2+[µγω 2 nr γ sin(πγ)]2 (2.17) In some cases of vibration analysis, K2n could be neglected in comparison with K1n (Kempfle et al., 2002; Rossikhin and Shitikova, 1997). In the case of aviscoelastic integer orderKelvin-Voigtmaterialmodel, thegoverning equation of the analysed beam has the following form (Kaliski, 1966) EJ (∂4y(x,t) ∂x4 +µ ∂5y(x,t) ∂x4∂t ) +Aρ ∂2y(x,t) ∂t2 =Fδ(x−g(t)) (2.18) The solution to the equation (2.18) could be obtained using a similar approach as in the case of the equation with fractional damping. For each n-th mode, the response could be calculated as follows (Kaliski, 1966) ξn(t)=    f0 ωhn t∫ 0 e−hn(t−τ) sin(ωhn(t− τ))sin nπg(τ) l dτ for hn <ωn f0 ω̃hn t∫ 0 e−hn(t−τ) sinh(ω̃hn(t− τ))sin nπg(τ) l dτ for h>ωn (2.19) where hn = 1 2 µω2n ωhn = √ ω2n−h 2 n ω̃hn = √ h2n−ω 2 n It should be noted that when γ = 1 (integer order derivative), the component K2n of frac- tional Green’s function (Eqs. (2.15) and (2.17)) vanishes, and it could be demonstrated that equations with a fractional and integer order derivative (Eq. (2.19)) are equivalent in the case of a subcritically damped system, i.e. ωn >hn. 1438 J. Freundlich 3. Calculation results and discussion In order to demonstrate the feasibility of the described above method, exemplary calculations of a beam subjected to a moving load have been performed. Responses of the beam with the fractional viscoelastic model are computed using equations (2.13)-(2.17) whereas the responses of the beamwith the integer order derivative model are computed using equation (2.19). The dimensionless dynamic deflection of the beam y/y0 versus the dimensionless time para- meter s has been computed. The variable y0 denotes the static deflection at the beammid-span and can be calculated from the equation y0 = F0l 3 48EI (3.1) The dimensionless time parameter s is defined as s= vt l or s= εt2 2l (3.2) in the case of a constant force velocity or constant acceleration, respectively. The calculations are accomplished for several values of the order of fractional derivative γ, dimensionless velocity v/vcr and acceleration ε/εcr. Variables vcr and εcr denote the critical velocity and acceleration, respectively. The critical velocity corresponds to the first natural vibration frequency of the beam and is defined as (Fryba, 1972) vcr = π l √ EJ ρA (3.3) whereas the critical acceleration is defined as the acceleration at which the force at the end of the beam reaches the critical velocity and is defined as εcr = π2 l3 EJ ρA (3.4) Firstly, the dynamic response of the beam subjected to a moving force at a constant velocity calculatedwith the help of presented above equations are compared to the results obtained using the corresponding formula given byFryba (1972). These comparative calculations are performed assuming a very light damping i.e. µγ = 1 · 10 −5 sγ and γ = 1 (integer order derivative). The calculations are performed for v/vcr = 0.5 (Eq. (3.1)). The calculation results obtained using the both equations are virtually identical. Minimal differences are found in the results, which is probably caused by differences in the modelling of damping in the system. Namely, in the equations given by Fryba (1972) it is assumed that the damping coefficient is independent of vibration modes and frequency in contrast to the damping model utilized in this work (Eqs. (2.3) and (2.13)). Next, in order to illustrate the application of the presented above procedure, computational examples for the response of the beam loaded with a moving force are perfor- med. The calculations are carried out for the beam of length 20m, mass density 7600kg/m3, cross-section area 2 ·10−3m2, cross-section moment of inertia 3.953 ·10−6m4, Young’s modulus 2.1 ·105MPa, coefficient µγ =3 ·10 −2 sγ. The force value is assumed F =100N. From equation (2.16), it follows that the values of damping coefficients and vibration ampli- tudes depend on the roots of the characteristic equation and, therefore, knowledge of the roots is required for further calculations. The roots are computed for the varying fractional derivative order γ. The dependence of the roots of the characteristic equation on the fractional derivative order is shown in Fig. 2 (for the first two vibrationmodes). The real part of the root determines the damping coefficient whereas the imaginary part determines the vibration frequency (Eq. Dynamic response of a simply supported viscoelastic beam... 1439 Fig. 2. Relationship between roots of the characteristic equation and the order of fractional derivative: (a) real part, (b) imaginary part (2.16)). Performed calculations reveal that the value of the damping coefficient increases with the order of fractional derivative (Fig. 2a) but the vibration frequency is practically constant (Fig. 2b). Moreover, the ratio of the real to the imaginarypartof the roots of the characteristic equation increases with the increaasing order of fractional derivative (Fig. 3a), whereas the amplitudeαn (Eq. (2.16)) is practically independent of the order of fractional derivative (Fig. 3b). It could be expected that the dynamic deflection of the beam decreases with the increasing order of fractional derivative. Fig. 3. (a) Ratio of the real to the imaginary part of the roots of the characteristic equation, (b) relationship between the amplitude α n and the order of fractional derivative The responses of the beam subjected to a force moving at a constant velocity are presented in Figs. 4-6. In Fig. 4, the dimensionless dynamic deflection of the beam at the point under the travelling force at various values of the derivative order γ is presented whereas the dynamic deflection at the mid-span of the beam is shown in Fig. 5. These results reveal that deflection of the beam under a moving force decreases with the increasing order of fractional derivati- ve γ (Fig. 4). The maximum deflection of the beam is greater for the dimensionless velocity v/vcr = 0.5 than v/vcr = 1.0. In the case v/vcr = 0.5, the dynamic deflection at the mid-span of the beam decreases with the increasing order of fractional derivative when s < 0.775 while increases when 0.775