Investigating Fractional Damping Effect on Euler–Bernoulli Beam Subjected to a Moving Load

. In this work, the dynamic response of Euler–Bernoulli beams of four diferent boundary conditions with fractional order internal damping under a traversing moving load is investigated. Te load is assumed to be moving with diferent values of constant velocity. A proposed approach to obtain the closed-form solution of the problem based on Green’s functions combined with a decomposition technique in the Laplace transform domain is introduced. Several cases are studied and compared to the literature; for instance, if simply supported beam is considered, the following three cases are to be explored: the case of elastic (or undamped) beam, the damped (or viscously damped) beam, and fnally the fractionally damped beam modeled by the fractional Kelvin–Voigt model. Te efects to the beam dynamic response induced by magnitude of moving load velocity, damping ratio, and fractional damping order are explored. Te results expressed sufcient agreement with similar problems found in literature and evidenced that the dynamic response of beams is signifcantly afected by varying the fractional order of beam damping as well as the moving load velocity. Accordingly, using fractionally damped materials exhibits better realistic behavior of beams and intermediate between elastic and viscous beam behaviors.


Introduction
Te subject of dynamic response under diverse types of loading for beams of viscoelastic materials is considered of utmost importance. Many applications of beams are found in bridges, railways, aircraft, and vehicles. As a result, various progressive studies have been performed to address the damping conduct of such materials (Rossikhin et al. [1], Paunović et al. [2], and Klanner et al. [3]).
Te prevalent use of the Euler-Bernoulli beam model is attributed to its simplicity and ability to give reasonable estimations for various engineering issues. Utilizing the Euler-Bernoulli's equation, Hilal and Zibdeh [4] presented a paper contributory to the major issue concerning the vibrations of an elastic homogenous isotropic beam excited by moving loads under generic boundary conditions. Te beam was subjected to moving force with constant velocity and according to the beam's response; closed-form solutions were acquired. In this study, Sumelka et al. [5] reformulated the classical Euler-Bernoulli theory utilizing fractional calculus. Such generalization was called fractional Euler-Bernoulli beams, and it resulted in nonlocal spatial description. After two years, Blaszczyk [6] derived the Euler-Bernoulli beam equation by using the variational approach that leads the equation to contain left and right fractional Caputo derivatives simultaneously. Ten, it was transformed into an integral equation to be solved analytically and numerically.
On the other hand, Fractional derivative models have proven to be of signifcant value for the precise characterization of the damping behavior of various types of materials (Bagley and Torvik [7]). Leibnitz and De L'Hospital were the frst to have introduced the concept of fractional derivatives in the 16th century (Podlubny [8]). When compared to diferent mathematical models including the Maxwell model, Kelvin-Voigt model, or other models instituted by Flügge [9], Pipkin [10], and Christensen [11], the fractional derivative model was found to coincide with the experimental data very well utilizing fewer parameters according to Di Paola et al. [12]. Tus, studying the damping properties of materials using these techniques is becoming essential and proven to be efcacious.
Adolfsson [13] used the fractional order model in investigations involving viscoelastic materials to obtain the large deformations of investigated beams. Studies on the dynamic response of beams resting on fractional viscoelastic foundation are also found in the literature. Te efect of the moving load velocity on such beams was presented by Ouzizi et al. [14], while the same authors investigated the efect of inducing multiple harmonic moving loads with a constant speed on the beam dynamical response (Ouzizi et al. [15]). More recently, in the study by Praharaj and Datta [16], the same beam response was investigated again when subjected to a moving load. In this study, the fractional order derivative-based Kelvin-Voigt model was used to describe the rheological properties of the viscoelastic foundation, while the Reimann-Liouville fractional derivative model was applied for a fractional derivative order, and the resulted fractional diferential equation of motion was solved by applying the modal superposition method and triangular strip matrix approach. In all of these studies, the efects of diferent system parameters on the beam response were evaluated.
Te solution of a fractionally damped beam equation is examined by Jena et al. [17] by applying the homotopy analysis method to calculate the dynamic response. Te unit step and unit impulse functions are deliberated for this analysis, and the acquired results for diferent orders of the fractional derivatives are compared well with others and achieved by the Adomian decomposition method, which was used in 2007 by Liang and Tang [18] to derive the analytical solution of viscoelastic that continues fractionally damped beam considering step and impulse function responses.
In 2016, a study involving a nonproportionally damped Euler-Bernoulli beam, which was subjected to moving loads, determined an analytical solution to assess the dynamic properties of such systems (Svedholm et al. [19]). In this investigation, it was assumed that the system's dynamic behavior can be assessed by superposition considering that adequate orthogonality conditions were derived and a closed-form solution for the dynamic behavior for a given eigenvalue was proposed. In the same year, Freundlich adopted the modal superposition approach to obtain the solution for force vibrations of a fractionally damped beam, while a convolution integral of fractional forcing function and Green's function were utilized to acquire the beam's behavior. Ten, an assessment of the beam's dynamic response under diferent fractional derivative orders and the moving load's velocities was performed. Green's function model was also used back in 1997 by Foda and Abduljabbar to investigate the efects of various parameters on the dynamic behavior of a plainly propped Euler-Bernoulli beam submitted under a transverse moving load. Te process was then proven to be simple and efective.
In this study, Labȩedzki et al. [20] modeled both stifness and damping terms in the Euler-Bernoulli beam equations using fractional derivatives for the Cantilever case this time. Te equations were formulated and solved for beams with and without tip mass. Te other two boundary conditions (namely, fxed-fxed and Cantilever-fxed) were considered by Blaszczyk et al. [21] to be analyzed by the fractional Euler-Bernoulli beam equation, where the diferential equation of motion was converted into an integral one considering the mentioned boundary conditions and exact solution was obtained, which contained a composition of the left and right Reimann-Liouville integrals. Te study presented three solutions for a constant, power, and trigonometric functions. On the other hand, Abro et al. [22] presented an analytic study of a simply supported beam case based on the modern fractional approaches utilizing Caputo-Fabrizio and Atanagna-Baleanu fractional diferential operators. Te equation of motion was fractionalized to investigate the efects of principal parametric resonances, and Laplace and Fourier sine transforms were invoked for investigating the exact solution.
Considering that the fractional derivative provides reliable models of fractionally damped structures and that few contributions have been introduced in this feld so far, the need is still there to realize more generic closed-form solutions and investigate the impact of diferent parametric changes on beams' response under diferent loading scenarios. Within this context, this study introduces a novel formulation of the closed-form solution of the Euler-Bernoulli beam with internal damping expressed by the fractional derivative, when traversed by a concentrated load moving with constant velocity, while considering the following four beam boundary conditions (BCs): pinned-pinned (PP), fxed-fxed (FF), fxedpinned (FP), and fxed-free or Cantilever (FC). In the introduced approach, the governing equation of the Euler-Bernoulli beam described using the fractional Kelvin-Voigt model is written as a fourth order partial diferential equation and Caputo's defnition is employed for the fractional derivative. Te orthogonality conditions are introduced to convert the partial diferential equation into a second order ordinary diferential equation, while at this stage, Laplace transform is exercised using the decomposition technique introduced in the study of Abu-Alshaikh et al. [23], combined with Green's theorem to obtain the closed-form solution of the problem.
For verifcation purposes, the obtained closed-form solution of the fractionally damped beam is utilized to obtain solutions for the beam without damping, and with integer order damping cases, these solutions are then compared with those available in the literature.
To conclude, the dynamic response of the fractionally damped Euler-Bernoulli beam having four diferent beam BCs (namely, pinned-pinned (PP), fxed-fxed (FF), fxedpinned (FP), and fxed-free or Cantilever (FC)) is investigated in this paper. Te analysis is based on Green's functions' approach combined with a decomposition technique in the Laplace transform domain; the closed-form solutions are analytically generated, and the results are discussed through the selected numerical examples. Furthermore, the impacts bestowed upon the beams' response by the moving load velocity, damping ratio, and the order of fractional derivatives are illustrated. It is worth mentioning that all the obtained solutions are truncated using the mathematical software Maple (Maplesoft [24]). Te resulted fractional derivative damping models may allow researchers to choose more suitable models to precisely ft experimental ones.

Basic Concepts and Formulations
Te fundamental concepts of a beam model with fractional damping are introduced in this section. Primarily, detailed presentation of the basic formulation for the fractional Euler Bernoulli beam is shown using the fractional Kelvin-Voigt model. Tis is followed by a step-by-step solution of the equation of motion for diferent case studies deploying relations for internal damping of viscoelastic materials and natural frequencies of vibration modes for elastic beams, the generalized Delta function, fractional calculus, Green's functions, and some other special functions' concepts.
In order to study the beam with the fractional model of internal damping, the elastic behavior of the material should be studied with the aid of the Kelvin-Voigt model for creep, (Christensen [11]), which reveals the relation between the stress σ(t) and strain ε(t) in the elastic region, and since the strain rate, _ ε(t), is appearing in that model, it is possible to implement the theory of the fractional order derivatives to write the stress variations with time for any value of the fractional derivative, α, and its value in our work is 0 < α ≤ 1. According to the fractional order Kelvin-Voigt model, which was proposed by Shermergor [25] and adopted in several following works, for instance, Rossikhin and Shitikova [26], and with further manipulations related to the derivation of this model, as in the study by Nasir [27], it is convenient to write the following: where E is the modulus of elasticity, μ is the viscous damping coefcient of the beam, C α is the characteristic internal damping coefcient of the material, and Γ (.) is Euler's Gamma function. Moreover, Caputo's defnition of the fractional derivative, Podlubny [8], is used to write the fractional derivative in the integral form, and it is going to be used in this work to solve the beam's fractional order equation of motion, as will be discussed during the analysis. After this brief discussion, as equation (1) illustrates, it can be found that one simple fractional model, the Kelvin-Voigt fractional model, can be applied for investigating the dynamics of viscoelastic materials, and this model has an intermediate behavior between elastic and viscous materials. For more discussion regarding this model and its applications, the interested reader is advised to visit (Bonfanti et al. [28]) and references therein.
Based on the aforementioned analysis, the governing equation of the Euler-Bernoulli beam described using the fractional Kelvin-Voigt model can be written as a fourth order partial diferential equation (Di Lorenzo et al. [29]) in the form as follows: where w(x, t) is the transverse displacement of the neutral beam axis at point with coordinate x and time t, m is the moving load, which has a constant velocity v, and the term δ(x − vt) represents the Dirac delta function, which has a unity magnitude when x � vt; otherwise, it is equal to zero. Moreover, the term D α 0 + (w(x, t)) � (z α w(x, t)/zt α ) represents the fractional derivative of order α. Tis term is commonly introduced in the literature using various defnitions; one of the most well-known defnitions is the Caputo fractional derivative, Podlubny [8], which can be written as follows: Te model of equation (2) is a general representation of viscoelastic Euler-Bernoulli beams that can have various BCs, and this model is used to solve the problem under consideration for various beam end conditions: the pinnedpinned, fxed-fxed, pinned-fxed, and fxed-free cases, and diferent values of α between zero and unity are to be considered. In obtaining the solution for equation (2), the density ρ, modulus of elasticity E, moment of inertia I, and the characteristic coefcient of the material C α are all assumed to be constants.
For any vibration mode, φ n (x) is considered, while Y n (t) is the vibration mode of the beam, the modal form of w(x, t) can be written in the form of infnite series as follows: In equation (4), the subscript n denotes the mode number investigated, and the accuracy of the results depends on the total number of modes truncated using the series; Shock and Vibration 3 only the frst mode of vibration is truncated in this work for which acceptable levels of accuracy are obtained. Tis is going to be discussed in detail during this analysis. Nevertheless, the theoretical procedure implemented enables the researcher, if desired, to efectively investigate a greater number of modes. Te undamped free vibration mode φ n (x) is written for general beam boundary condition as follows: where A n , B n , and C n are constants evaluated from the BCs and c n is the root of the frequency equation, and all are obtained as listed in Table 1.
When substituting equations (4) into (2) and considering that C α � E(μ/E) α , we get the following: Multiplying each term by φ m (x) and integrating over the entire length of the beam from 0 to L gives the following: Now, it is possible to apply the orthogonality conditions, which imply the following important results: Using trigonometric and hyperbolic integration relations, it is simply shown that equation (9) is applicable for any combination of BCs (Rao [30]). Tus, in terms of equations (8) and (9), (7) becomes Let P 0 � mg, then the diferential equation of the n th mode of vibration of the generalized displacement or modal response of the beam may be rearranged and written as follows: Ten, equation (12) can be solved using a decomposition method similar to that used in the study of Abu-Mallouh et al. [31] for a linear beam equation of motion; nevertheless, if the equation of motion is nonlinear, the same decomposition is primarily used and then another decomposition is applied to resolve the nonlinear terms; this process is normally called the Adomian decomposition method (ADM) for nonlinear systems, refer to Kwong et al. [32] for more details regarding the ADM and its applications. From equation (12), the solution of the vibration mode of the beam, Y n (t), is assumed to have a series form as follows: So, equation (12), in view of equation (13), is written as follows: Substituting the n th natural frequency ω n and the damping ratio ζ n as follows: Tus, equation (14) becomes and by rearranging equation (16), we get the following: where α n � sin(c n l) + sin(c n l)/ cos(c n l) + cos h(c n l) Shock and Vibration 5 Applying the decomposition form illustrated in equation (13), equation (17) can be divided into a system of two recursive equations as follows: It is noteworthy that in equation (18), the subscript n denotes the vibration mode to be studied, while in equation (19), the superscript J is used to indicate the number of terms considered in the decomposition. In the following equations, for enhanced readability, the summation signs in equations (18) and (19) are dropped (same for the following equations); nevertheless, when writing the fnal solutions, the summation terms are to be reconsidered. Now, by defning Φ n (s) to be the Laplace transform of the vibration mode function, as shown in equation (5), then we get the following: where Ω n � c n v. Taking the Laplace transform for equation (18), while assuming homogenous conditions, yields to the following: Ten, the initial part of the solution in the Laplace domain is simplifed by rearranging equation (21) in terms of y 0 n (s) and then applying the direct Laplace inverse to get the initial solution in the time domain denoted by Y 0 n (t), as follows: Taking the Laplace transform for equation (19), while assuming homogenous conditions and considering that J ≥ 1, yields to the following: Rearranging equation (23), the recursive formula for the solution is obtained as follows: So, using equation (24) with the aid of the value y 0 n (s) from equation (21) and applying math induction, the general form of y J n (s) for J ≥ 1 is written as follows: Ten, the total series solution, in the Laplace domain, is written as follows: Tis can be rearranged as follows: Comparing the series in equation (27) with the general form of the geometric series, that is, In this, Υ � P 0 Φ n (s)/χ n (s 2 + ω 2 n ) and R � −2ω n ζ n s α / s 2 + ω 2 n . Realizing that the quantities ω n , ζ n are positive values, then if 2ω n ζ n s α /s 2 + ω 2 n < 1, the summation presented in equation (27) according to equation (28) becomes the following: As a result, to fnd the solution in the time domain, it is desired to fnd the Laplace inverse for equation (29), and for this purpose, a number of cases will be studied hereafter.

Elastic Beam.
Beams behave elastically when the damping ratio equals zero; therefore, the closed-form solution of the elastic beam can be written as follows: 6 Shock and Vibration n ω 2 n + Ω 2 n Ω n sin ω n t − ω n sin Ω n t + ω n A n ω 2 n + Ω 2 n cos ω n t − cos Ω n t + B n ω 2 n − Ω 2 n Ω n sin ω n t − ω n sin h Ω n t To verify the validity of this solution, the simply supported beam case is considered, and the coefcients of the vibration mode in equation (30) are all set to zero, then the solution becomes the following: Tis equation represents exactly the same result obtained by other formulations made in the literature (Rao [30] and Abu-Malloh et al. [31]). Finally, by substituting equation (30), with equation (4), the transverse displacement of the system w(x, t) at point x and time t for the elastic beam can be written in the following form:

Integer Order Damped
Beam. Considering the viscous beam described by the Kelvin-Voigt model while being affected by a moving load with constant velocity, then the generalized displacement of this special case can be obtained as where Y n1 (t) � P 0 χ n Z 2 1 Z 1 Ω n 2ω 2 n ζ 2 n + Ω 2 n − ω 2 n sin h Z 1 t cos h ω n ζ n t − sin h ω n ζ n t + 2ω n ζ n Ω n cos h Z 1 t cos h ω n ζ n t − sin h ω n ζ n t − cos Ω n t − Z 3 sin Ω n t Shock and Vibration 7 in which Finally, by substituting equation (33) in equation (4), then w(x, t) for the beam with integer order damping can be written in the following form: Y n (t) sin c n x + A n cos c n x + B n sin h c n x + C n cos h c n x . (36)

Fractional Kelvin-Voigt Model.
In this section, the fractional Kelvin-Voigt modeled beam is studied; the beam is assumed to be fractional with order of α, where 0 ≤ α ≤ 1, referring to equation (29), and the recursive formula for this beam in the Laplace domain can be written as follows: For the frst term, y n0 (s), the direct Laplace inverse is used to fnd the solution, as it was previously achieved in equation (22). Ten, Laplace inversion for y n1 (s) term can be obtained by applying the convolution theorem as follows: y n1 (s) � 2ω n ζ n P 0 s α Φ n (s) χ n s 2 + ω 2 n s 2 + 2ω n ζ n s α + ω 2 n � 2ω n ζ n P 0 Φ n (s) χ n s 2 + ω 2 n ⎡ ⎢ ⎣ ⎤ ⎥ ⎦ s α s 2 + 2ω n ζ n s α + ω 2 n � f(s) g(s). (38) For f (s), Laplace inverse is directly obtained as follows: For g(s), the following fractional Green's function of the third order is used to obtain the Laplace inverse (Podlubny [8]) as follows: Ten, g(s) can be rewritten as follows: Te coefcients are as follows: Hence, G(t) becomes 8 Shock and Vibration where E (k) 2,2−α+αk (−ω 2 n t 2 ) is the k th derivative of the Mittag-Lefer function, which is defned as follows: So, equation (40) becomes the following: Now, applying the convolution theorem for the functions presented in equations (39) and (45) leads to Te solution of the integration in equation (46) is found analytically by rearranging equation (45) as follows: Defning the following parameters, we have Hence, G (t) becomes the following: By computing integration in equation (46) term by term, the value of Y n1 (t) is found as a summation of the following four equations: Hω n t H+1 LommelS1 1/2H, 3/2ω n t Tus, where the hypogeom(.) and LommelS1(.) are the generalized hypergeometric and generalized Lommel special functions, respectively, and then Y n (t) is obtained by adding the two solutions in equations (22) and (52). Finally, w(x, t) for the fractionally damped beam can be written in the following form: w(x, t) � (53) All the solutions of the previously presented formulations are truncated and fulflled using the mathematical software package, Maple.

Results and Discussion
In this section, a numerical verifcation will initially be presented to verify the obtained formulation in the previous section, and then the application of these equations is extended to plot the generalized defection response of the beam to a moving load for the four BCs considered. Moreover, the proposed results are to be compared to those corresponding results present in the literature obtained using diferent approaches. Finally, the efects of the load normalized velocity, damping ratio, and the order of the fractional derivative on the beam's response are presented and discussed.

Numerical Verifcation of the Formulations.
In order to perform the numerical verifcation of the solutions of equations for the frst mode shape (n � 1) and after considering that the nondimensional length of the beam is unity, some arbitrary set of parameters are taken as follows: Tese parameters are to be substituted in the corresponding equations to obtain the response of the beam in terms of maximum defection. It is important to mention that, for most results shown in this section, a nondimensional analysis is practiced, for which the units of the above parameters are eliminated, and this practice is advantageous to generalize the results for any desired set of data.
Firstly, the response of the elastic beam case in terms of maximum defection at midpoint of the PP beam due to the moving load is calculated by the closed-form solution, as shown in equation (30). Te resulted values of this solution are considered as control. Ten, the beam with integer order damping case equations are used to calculate the response of the elastic beam, which can be achieved by setting the value of ζ to be approaching zero. Systematically, the elastic beam response is recalculated by the fractional Kelvin-Voigt beam equations. In this case, both αand ζ are set to be approaching zero. Both results from integer order damped and fractional order Voigt equations are compared to the control solution as listed for each term in Table 2.
It can be seen from Table 2 that, for the simple case of elastic beam, all equations obtained within this contribution provided almost the same solutions with an error approaching zero percent. It is worth noting that results from the frst mode of vibration are used in the whole paper because it is proved to be the most efective one in such beam cases.
Te same procedure of verifcation is performed this time by solving integer order damped PP beam case using equation (33) by setting ζ � 0.25 and then comparing the resulted defection at beam midspan to those obtained using the fractional Voigt beam equations when α � 1 and ζ � 0.25, which in this case represents a beam with integer order damping. Table 3 shows that the solution of defection obtained from the integer order damping solution fts exactly into the solution obtained from the fractional beam equations, which implies that both elastic and integer order damping beam cases can be solved by the proposed set of fractional beam equations.

Efect of the Moving Load Velocity.
In order to generalize the representation of the solution of the beam, dimensionless parameters are introduced, namely, dimensionless speed (β) of the moving load, which is the ratio of the load speed (v) to the load critical speed (C cr ), and the dimensionless time (t), which equals to tv/L. Here, C cr � ω 1 L/π as mentioned in the study by Hilal and Zibdeh [4]. Te normalized beam defection then can be defned as follows: where x max and (STD) max are location and value of the maximum static defection of the beam, respectively, and the normalized position x is defned as the position along the beam (x) divided by the beam length (L), as x � x/L. In Figure 1, the normalized elastic beam defections versus the dimensionless time are plotted for PP, FF, FF, and FC beams, considering fve diferent load velocities: β � 0.05,0.2,0.5,0.8, and 1. Te curves of these fgures conform well with those obtained in the literature considering similar parameters (Foda and Abduljabbar [33] and Fryba [34]). Moreover, from Figure 1, it can be clearly seen that at lower velocity, the response reaches its peak at a certain value of t and then decreases to zero at t � 1. While at higher velocities, the peak is shifted to the right, and it is still high at t � 1. Moreover, the response of the beam increases by increasing the velocity, until it reaches its peak at β � 0.5, and then the maximum response values start to decrease as β increases. On the other hand, at very low velocities, the response behavior is reciprocating since the moving load is not able to deform the beam efectively, while at higher velocities, the response behavior is more obviously afected by the moving load. In order to compare the response of the four types of beams under the moving load, β � 0.5 is chosen   as a speed parameter since the response at this velocity is relatively high compared with other velocities. Table 4 lists the maximum normalized defection of the four types of beams at β � 0.5, as shown in Figures 1 and 2, which show a comparison for the actual (dimensional) beams' response along the normalized beam length. From Table 4 and Figure 2, the maximum beam defection occurs earlier, in terms of time and distance, in the PP beam than the other cases.

Efect of Damping Ratio.
In general, when increasing the damping ratio, the response of the beam is normally decreased; to check if this applies on the results of the proposed formulations and to compare with the previously obtained results by other approaches, Figure 3 is plotted to show the response of the beam with integer order damping under varying damping ratios (ζ � 0, 0.1, 0.2, and 0.3). Te curves of these fgures ft exactly into the results obtained in the literature by Hilal and Zibdeh [4]. Te efect of the damping ratio is clear in this fgure; as for ζ � 0.0, the beam demonstrated exact elastic response, and when increasing ζ, responses get closer to each other, while at t � 1, all the responses are almost zero. Figure 4 shows the dynamic response of the four beams when β � 0.5 and ζ � 0,0.1,0.2, and 0.3. Te curves of this fgure confrm that when ζ � 0, beams respond elastically; however, when increasing the damping ratio, the peak defection decreases and gets shifted to a higher value of t. Also, the defection of the beam in the fgures decreased when the damping ratio increased, and this result agrees with the concept of damping since the damping reduces the vibration efect of the beam. Table 5 lists the percentage reduction of the peak normalized defection for each type of beam under study, in which the PP case exhibited more reduction percent in its normalized defection than the other cases.
An attempt to study the efect of the moving load velocity on the beam response at diferent damping ratios is made by plotting Figure 5. Tis fgure, in particular Figure 5(a), exhibits that the same trend appeared in similar fgures found in the literature [4], while Figure 6 compares the response of diferent end condition beams when they are damped with integer order at β � 0.5 and ζ � 0.1.

Efect of the Fractional Derivative.
In this subsection, the efect of the beam's fractional derivative order to the beam's response is investigated, and this efect is not widely studied in the literature, especially for general beam BCs. Nevertheless, for the simply supported beams' case, some previous works [35] are found and compared to the results of the present work. Figures 7 and 8 are plotted to demonstrate the efect of a wide range of varied velocities of moving load on the midspan responses of the PP and FF beam, respectively, when fractional order α � 0.5, while in Figure 9, the fractional derivative α is varying and the responses of the beam to a moving load of β � 0.5 are plotted. Te damping ratio for all of these fgures is taken as ζ � 0.1. When comparing the curves of Figure 7 to the curves of the available cases in the      literature (β � 0.25, 0.5, 0.75, and 1 only), it can be said that this fgure presents a solution that conforms to others presented previously (Freundlich [35]), while the curves of Figure 9 showed equal peaks' response to similar curves in the literature; however, diferences appeared due to the numerical solution accuracy and to the efect of the damping ratio considered. Figure 9 implies that, by decreasing the order of the derivative, the response of the beam increases; nevertheless, the peak of the curve is not shifted, so the beam still reaches its peak value at the same time. Moreover, the noticeable efects of the fractional derivative on these curves are around the normalized peaks' region and at the end of the normalized time region rather than the rest of the curve regions. If the curves of Figures 7 and 8 are compared with the beams damped with the integer order, clear diferences in the peak values are noticed; hence, decreasing the order of the fractional derivative yields to increasing the normalized defection of the beams. In Figure 10, a comparison among the responses of fractionally damped beams of the four end conditions is presented at β � 0.5, ζ � 0.1, and α � 0.5. Table 6 supported by Figure 11 presents the maximum dimensional defection of the elastic, integer order damped, and fractionally damped beams at diferent BCs, while ζ � 0.1 and β � 0.5. Te data in Table 6 show that for beam types studied, when increasing the order of the fractional derivative, defection decreases; moreover, the response of the fractionally damped beams is always lower than the elastic beam and greater than the damped beam with integer order; that is, the behavior of the fractionally damped beam is intermediate between the elastic and fully viscous beam.

. Conclusions and Recommendations
Tis study investigated the dynamic response of Euler-Bernoulli beams of four diferent BCs with fractional order internal damping traversed by a moving load. An approach is introduced to obtain the closed-form solution of the problem based on Green's functions combined with the decomposition technique. Comparisons with the previous studies available in the literature for verifcation purposes showed great matching, and based on these verifcations, further results were performed for cases that have not been available in the literature yet. Te moving load efects were   investigated for a range of diferent values of constant velocity. In addition, the efects of the beams' dynamic response caused by the order of fractional derivative internal damping were also explored.
Te main conclusion of this research is that the efects of using the fractional order derivative model containing fractional constitutive law gives some good consequences, despite the difculties of using the fractional model due to the complicated math techniques and programming eforts.
From the obtained results and fgures, the following conclusions can be also drawn: (a) Te amplitude of normalized defection increases with the moving load velocity increase to a certain value, and then it starts to decrease. (b) Te beam dynamic oscillations were more evident for lower values of moving load velocities. In particular, PP beam sufered clearer oscillations. (c) Positions of the beam's amplitude response shifts to the right when the load velocity increases. (d) Among the investigated load velocities, the maximum beam normalized response appeared at β � 0.5 for all BCs. (e) When damping ratio was decreased, beams sufered more oscillations and response amplitudes were shifted to the left. (f ) Te largest amplitude of normalized defection appeared for PP case and the lowest appeared for FC, while the largest amplitude of dimensional response appeared in the FC case and the lowest was for FF. (g) Te fractional order afected most of the amplitude response of FF and the least that of FC, while for PP and PF, efects were clearer at the time end. (h) Beam defection at α � 0.5 is found to be less than the average defection of those when α � 0.0 (undamped case) and α � 1.0 (one integer damping case), which implies that the defection does not change linearly with the order of the fractional derivative α.
For any certain engineering application, further future research can be implemented, in order to investigate the optimal fractional order derivative (α) that yields to the best desired beam defection. Moreover, it has been shown that, for the frst natural frequency, the proposed analytical model results in excellent beam defection results. Nevertheless, for other beam's natural frequencies, validations of the proposed analytical model could be performed. Moreover, the impact of various mass ratios between the moving load and the total weight of the beam could be held in future studies. Additionally, the problem can be improved for moving vehicles, trains, or multiloads or masses, and the efect of the functionally graded material on the dynamic response may also be studied. Furthermore, the efect of large deformations on the beam dynamic defection can be studied, where a nonlinear term is to be added to the beam governing equation, and in this particular case, approximate numerical solutions are expected to be generated.

Data Availability
Te data that support the fndings of the study are available from the corresponding author upon request.

Conflicts of Interest
Te authors declare that they have no conficts of interest.