Transient vibrations of a fractional Zener viscoelastic cantilever beam with a tip mass

The presented work concerns the kinematically excited transient vibrations of a cantilever beam with a mass element fixed to its free end. The Euler–Bernoulli beam theory and the fractional Zener model of the beam material are assumed. A fractional Caputo derivative is used to formulate a viscoelastic material law. A characteristic equation, modal frequencies, eigenfunction and orthogonality conditions are achieved for the beam considered. The equations of motion of the system are solved numerically. A numerical solution of a multi-term fractional differential equation is obtained by means of a conversion to a mixed system of ordinary and fractional differential equations, each of the order of 0<γ≤1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$0 < \gamma \le 1$$\end{document}. The transient time histories of the beam vibrations during the passage through resonance are calculated. A comparison between the beam responses obtained with a fractional and an integer viscoelastic material model is presented. The calculations performed reveal that use of the fractional damping affects on the time histories of the system. The calculated beam responses show that for some values of the order of the fractional derivative γ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\gamma$$\end{document}, the amplitudes occurring in the area of the second resonance are greater than those obtained in the area of the first resonance, which does not occur in the case of the integer order of the fractional derivative. Moreover, an evaluation is made of the difference between the results obtained for the calculations using the fractional Zener model and the fractional Kelvin model. It is shown that for some physical beam parameters, the calculation results obtained using both models are virtually the same for both models, which means that the the simpler, fractional Kelvin–Voigt material can be used instead of the fractional Zener material model. This simplifies the solution and decreases the time needed to make the numerical calculations.


Introduction
Its well known that cantilever beams are widely used to model various structural elements such as masts, offshore structures, military airplane wings, accelerometers, Stockbridge dampers, energy harvesters and high buildings [1][2][3][4][5]. A number of authors have studied the dynamics of a cantilever beam with a mass element fixed to its free end [6][7][8][9][10][11][12][13]. The majority of these studies consider only the beam dynamics without damping or with a viscoelastic material model, where an integer order Kelvin-Voigt material model is employed. Furthermore, these works seldom examine the forced transient vibrations of the analyzed system. Machine components often work in an oscillating environment, which can cause undesirable vibrations of these elements. Moreover, machines frequently operate in an area above their first or subsequent resonance and therefore have to pass through one or more resonances when starting or stopping them. Therefore, it is very important to determine the maximum amplitudes of the vibrations and stresses that occur in a system while it passes through a resonance. Many authors have studied this issue for various system configurations, but many such studies leave out some important considerations, such as the influence of damping on system dynamics and orthogonality conditions (see e.g. [3,[6][7][8]). It is well known that the viscoelastic properties of materials can significantly affect a system's dynamic behavior. An important issue in modeling the dynamics of a structure is making a correct determination of its energy dissipation properties, for the system's dynamic behavior is strongly dependent on its damping properties, which also depend on the viscoelastic properties of the material. Thus, a dynamic analysis should involve an appropriate model of the viscoelastic material. Many structural materials show weak frequency dependence for their damping properties over a broad frequency range. This weak frequency dependence of damping requires a complex description of the viscoelastic material using differential equations of the integer order [14][15][16][17]. In recent decades, fractional differential equations have been used increasingly frequently to describe many physical phenomena, including the properties of viscoelastic materials [18][19][20]. The fractional derivatives are used to model both spatial (e.g. [21][22][23][24]) and time dependent physical phenomena [18,20,25]. By means of fractional derivatives it is possible to model the damping properties of viscoelastic materials, especially those that show a weak dependence of the damping properties on the frequency [14,17,26]. The constitutive equations describing viscoelastic materials whose damping properties are weakly dependent on the frequency, written with fractional differential equations, require fewer parameters than those written using integer order equations [15][16][17]. An overview of works and events concerning the fractional calculus is provided in a paper by Machado at al. [19], whereas a review of the historical development of fractional calculus in the mechanics of solids is presented by Roshikhin [27]. A survey of publications on the application of fractional calculus in the dynamics of solids is presented in a paper by Rossikhin and Shitikova [28].
Probably the first paper dealing with the transient dynamic response of a beam consisting of a viscoelastic material described by a fractional derivative model was published by Bagley and Torvik [28,29]. They analyzed a simply supported Euler-Bernoulli beam with a viscoelastic layer described by a fractional Kelvin-Voigt model. In their paper, the transient response of the beam to step loading using a continuum formulation and the finite element method was presented. Later, the dynamics of beams with viscoelastic properties described by the fractional Kelvin-Voigt model were studied in various works [28,[30][31][32][33]. It is the most commonly used fractional model of a viscoelastic material in modeling the dynamics of beams. The fractional Kelvin-Voigt model has some advantages, of which the most important is that it can be easily applied in the analysis of the dynamics of various structural elements. The use of the Kelvin-Voigt fractional viscoelastic material model in modeling the dynamics of structural elements allows us to obtain equations that can be relatively easily solved by known methods used in solving elastic problems. Unfortunately, this model also has some disadvantages. One of them is that it overestimates material damping for higher frequencies. Caputo and Mainardi [15] have show that the fractional Zener solid model allows for a wider agreement with the experimental data then fractional Kelvin-Voigt model. Similar conclusions have been presented by other authors, e.g Bagley [34] and Pritz [35]. Another disadvantage is that this model cannot describe the behavior of some types of viscoelastic auxetic materials possessing negative Poisson's ratios [36]. The fractional Zener solid model does not have the aforementioned disadvantages [36] and is therefore being used increasingly in modeling the dynamics of solids. The fractional Zener solid model has been employed in beam dynamics. A brief overview of the research using this material model to describe the beam dynamics is presented below.
Stankovic and Atanackovic analyzed the dynamics and stability of a beam with the fractional Zener material model. They studied the lateral vibrations of a compressed simply supported beam. They analyzed the existence of a solution in the space of Laplace hyperfunctions and performed a detailed stability analysis of the solution [37]. In the paper [30], aforementioned authors analyzed a similar beam system. They showed that the dynamics of the problem is governed by a single linear differential equation with a fractional derivative. Moreover, they obtained the form of the solution for the case of a periodic compressive force. Lewandowski and Baum [38] presented a dynamic analysis of a multilayered beam with viscoelastic layers described by the fractional Zener model. They applied the Euler-Bernoulli beam theory to elastic layers and the Timoshenko theory to the viscoelastic layers. They adopted the finite element method, the virtual work principle, and the Laplace transformation to derive of the equation of motion in the frequency domain. Additionally, the authors found that system parameter changes impacted greatly on the value of the natural frequency and the non-dimensional damping ratios of the beam. Lewandowski and Wielentejczyk [39] studied non-linear, steady state vibrations of beams, excited by harmonic forces. They applied the von Karman theory to take account of geometric nonlinearities. They used the finite element approach. The one-harmonic function in time was assumed to describe the steady-state vibration of the beams being studied. The authors utilized the harmonic balance method for a system of which the constitutive stress-strain relationship was given as a differential equation containing the fractional derivatives for both stress and strain. Furthermore, they examined the stability of the steady-state solution. Martin [40] presented quasi-static analysis, classical and fractional dynamic analysis of a simply supported viscoelastic beam subjected to uniformly distributed load. In this paper, the author studied the dynamic behavior of Euler-Bernoulli beam with its damping properties described in terms of fractional derivatives of an arbitrary order. The governing equation was solved with a mixed algorithm based on Galerkin's method for the spatial domain and the Laplace transform, Bessel functions theory and binomial series expansion for the time domain. The proposed approach was used to solve sample problems. The usefulness of the proposed approach has been demonstrated by calculations of beam deflections in time for various orders of the fractional derivative.
Analyzing the work described above, it can be noted that the use of the fractional Zener solid model significantly complicates the solution of the investigated problem as well as further numerical calculations. In addition, the use of the fractional Zener model in the description of the forced vibration of structural elements causes the appearance of a fractional derivative of the forcing force in the equations of motion [41,42], which significantly complicates the solution of the problem and the numerical calculations. Furthermore, the presence of a fractional derivative of the forcing force in the equations of motion can be expected to significantly increase the numerical computation time. Thus, although the fractional Zener model better describes the dynamic properties of materials over a wide frequency range than the Kelvin-Voigt model does, its application greatly complicates the solution of the investigated problem and the numerical calculations.
In the author's opinion, it should be examined whether, in calculating the dynamics of beams or other structural members, using the fractional Zener model yields significantly different results than using the fractional Kelvin-Voigt model. Moreover, the number of research works that have studied the transient forced vibrations of viscoelastic beams described by the fractional Zener model is rather limited. In particular, there is a lack of studies of the transient forced vibrations of cantilever beams with a tip mass element and using the fractional Zener material model.
Analysis of the transient vibration of the cantilever beam with tip element mass and with material properties described by the fractional Kelvin-Voigt model has been presented in a previous paper by the author [32]. However, in this work the equations of motion were obtained without using the Laplace transform. Moreover, a different method of computing the generalized coordinates was used. Namely, the beam response was calculated using a convolution integral of the fractional Green's and forcing functions.
To the author's knowledge, a transient dynamic analysis of a cantilever beam with a mass element fixed to its end has been performed only for beams with a fractional viscoelastic material described with the help of the fractional Kelvin-Voigt model. The objective of this study, therefore, is to analyze the transient vibrations of a cantilever beam with a mass element fixed to its free end and with a viscoelastic material characterized by the fractional Zener model. Furthermore, the purpose of this paper is to investigates whether the use of the fractional Kelvin-Voigt model can be sufficient when calculating structural member dynamics, that is, whether the calculation results obtained using the two models mentioned differ significantly. The presented analysis may be useful, for example, in modeling energy harvesting devices [4].

Problem formulation
In this work, a cantilever beam with a tip element of mass m p , moment of inertia I B and length l is studied. A homogeneous beam with a uniform cross-section and mass density is assumed. The mass center of the mass element coincides with the free end of the beam (Fig. 1). The Euler-Bernoulli theory is assumed, i.e. the rotary inertia and shear deformation are neglected. It is assumed that the beam can move only in the xzplane and that the gravitational force is perpendicular to this plane, so that the gravitational force does not affect the beam motion. The beam is excited by a base motion along the z axis (Fig. 1).The viscoelastic properties of the beam material are assumed to be characterized by a fractional Zener model, which is defined as follows [15,16,43] rðtÞ þ bD ðcÞ ðrðtÞÞ ¼ E ðtÞ þ l c D ðcÞ ððtÞÞ ð1Þ where rðtÞ and ðtÞ are the stress and strain histories, b and l c are time constants [43], E is the relaxed modulus, t is time and D ðcÞ ðÁÞ is the Caputo fractional derivative of the order c defined as [14,18,44] D ðcÞ ðf ðtÞÞ 1 Cðm À cÞ where Cðm À cÞ is the Euler gamma function [18], D ðnÞ f ðÁÞ ð Þ ¼ o n ot n ðÁÞ is the n-th derivative of a function f ðÁÞ with respect to time, m is a positive integer number satisfying the inequality m À 1\c\m, and t [ 0.
It is often assumed that the order of the fractional derivative is in the range 0\c 1 for many real materials [16,17], where c ¼ 1:0 corresponds to the integer order derivative [18,44].
Using the fractional Zener material model Eq. (1), the following equation for the forced transverse motion of a beam is obtained (see ''Appendix A'' section) where w(x, t) is the transversal displacement of the neutral beam axis ( Fig. 1), A is the area of the crosssection of the beam, J is the moment of inertia of the beam cross-section with respect to the neutral axis, q is the material mass density of the beam, q(x, t) is the transverse load per unit length acting on the beam, x is the longitudinal coordinate. The solution is sought on the assumption that the beam is at rest at the initial instant of time, and so the initial shear force and the bending moment are zero. The solution to Eq. (3) must satisfy the boundary conditions for the beam considered, i.e. the displacement follows from a base motion whereas the slope of the displacement curve must be zero at the fixed end, i.e. x ¼ 0, and so the boundary conditions are defined as below where w bs t ð Þ is a function of time, describing a base motion.
The boundary conditions of the beam at the free end with the mass element are obtained using the rate of change of the linear and angular momentum of the mass element relative to its mass center (see Appendix ''B'' section) [45,46]. Thus, where V(l, t) and M(l, t) are the shear force and bending moment, respectively. The relationship for a bending moment can be obtained by multiplying both sides of Eq. (1) by z and integrating over the cross-section area of the beam (e.g. [38]) Mðx; tÞ þ bD ðcÞ ðMðx; tÞÞ For a beam excited by a base motion, the transverse translation of the beam can be considered as the sum of the rigid body and the relative translation [47], viz where w rl is the relative translation of the neutral beam axis with respect to the clamped end of the beam. Substituting appropriate derivatives of Eq. (7) Similarly, by substituting derivatives of the expression Eq. (7) into Eq. (6), the bending moment is expressed as Mðx; tÞ þ bD ðcÞ ðMðx; tÞÞ Assuming zero initial conditions, and applying the Laplace transform [18,44] to the equation above, the expression for the bending moment in the Laplace where f ðnÞ ðÁÞ ¼ o n ox n f ðÁÞ is n-th derivative of function f ðÁÞ with respect to the spatial variable, s is the Laplace variable.
Using the relationships for bending moment Eq. (10) and applying the Laplace transform to the boundary conditions for the end with a tip mass element Eq. (5), the boundary conditions at x ¼ l can be formulated as follows Next, the equations above can be transformed as Applying the Laplace transform to the equation of motion Eq. (3), one obtains this equation can be rewritten as where where The solution to the problem formulated in equations (12) and (13) is sought as a convergent series of the eigenfunctions of the elastic beam [28,45,48] b where W n ðxÞ are the eigenfunctions of the corresponding elastic beam, b n n s ð Þ are time-dependent generalized coordinates in the Laplace domain.
The eigenfunctions W n ðxÞ may be found using a well-known procedure (see e.g [4,32,45]), namely by solving analogical equation to Eq. (15) for the elastic beam with homogeneous boundary conditions Eq. (12). Then, using separation of variable method, the following equation can be obtained The solution to Eq. (18) is sought in the form Utilizing the boundary conditions for the clamped beam end, it can be found that A ¼ ÀC and B ¼ ÀD.
Next, employing these results and the boundary conditions for the beam end with the element mass, the system of equations for the constants A and B is formulated as (20) is fulfilled if the determinant of the coefficients matrix of the system of equations equals zero. Thus, equating the determinant of Eq. (20) to zero, after many mathematical transformations, the characteristic equation of the system is obtained The solution of the characteristic equation (21) has a countable infinite set of roots. The p n root of the characteristic equation corresponds to the n-th natural undamped frequency x n of the elastic beam. The frequencies can be calculated using the well-known expression [45,46,49,50] x n ¼ k 2 where Having calculated the roots of the Eq. (21), the numbers k n and the natural frequencies x n can be calculated. Next, determining k n from Eq. (22) and substituting into Eq.(16), the equation for natural damped frequencies of the viscoelastic beam is obtained [49] x 2 The equation above may be rewritten as In a paper by Rossikhin and Shitikova [51], it is shown that Eq. (24) has no real negative root and lacks complex conjugate roots in the right half-plane of the complex plane. The real part of the n-th root determines the damping coefficient of the n-th mode, while the imaginary part determines the n-th damped natural frequency. Substituting the roots p n and the numbers k n into Eqs. (20) and (21), the eigenfunction of the analyzed beam can be written as W n ðxÞ ¼ A n ½sinðk n xÞ À sinhðk n xÞ À k n ½cosðk n xÞ À coshðk n xÞ f g where k n ¼ p n Á a Á sin p n ð Þ À sinh p n ð Þ ½ À cos p n ð Þ þ cosh p n ð Þ ½ p n Á a Á cos p n ð Þ À cosh p n ð Þ ½ þ sin p n ð Þ À sinh p n ð Þ ½ Functions W n ðxÞ must fulfill the orthogonality conditions [48]. These orthogonality conditions have the following form [4,32] or, in an alternative form where d mn is the Kronecker delta, and m l ¼ Aq The eigenvalues and eigenfunction can be computed using equations Eqs. (21) and (25). The knowledge of these values enable us to obtain a solution to the forced vibration. Substituting Eq. (17) into Eq. (13), the following expression is derived W n x ð Þn n s ð Þs 2 1 þ bs c ð Þ ¼ qðx; sÞ À qAs 2 w sb ðsÞ Multiplying both sides of Eq. (28) by the massnormalized eigenfunction W s ðxÞ and integrating over the length of the beam, one obtains Employing the orthogonality conditions (Eqs. (26)- Arranging and grouping the appropriate terms, after several arithmetical transformations, the following relationship is obtained In the next step, the boundary conditions (Eq. (12) It can be noticed that the expressions under the sum operator in Eq. (31) are the same as the left-hand sides of the equations for the boundary conditions (Eq. (32)), and so, after a number of arithmetical transformations, the following equation for generalized coordinates in the Laplace domain is derived Assuming zero initial conditions, the equation above corresponds to the following equation in the time domain [18,44] D ð2þcÞ n n t ð Þ ð Þþ The differential equation achieved for the generalized coordinates for the fractional Zener beam (Eq. (34)) is more complicated than the analogous equation obtained for the fractional Kelvin-Voigt beam [32]. It can be seen that the order of Eq. (34) is 2 þ c, meaning that the order is fractional and it is greater than 2, which greatly complicates the numerical solution of this equation. Moreover, the right side of the Eq. (34) contains fractional derivatives of the load and inertial force, which also complicates the numerical solution and significantly extends the numerical simulation time.
Generalized coordinates n n ðtÞ can be calculated numerically with the help of a method similar to the methods described in the paper by Edwards et al. [52] and in the book by Diethelm [53] (Chapter 8). Namely, the numerical solution of the linear multi-term fractional differential Eq. (34) can be obtained by converting this equation to a system of mixed ordinary and fractional differential equations, each of the order 0\c 1. Therefore, the Eq. (34) can be converted into a system of equations as shown below where Q n ðx; tÞ and P n ðtÞ have the same meaning as in Eq. (34). It can be noticed that The identity above simplifies the solution of the system of differential equations Eq. (35) Due to the presence of the fractional derivative in the equations of motion, the use of the method presented in the author's previous paper [32] was practically impossible, due to the considerable difficulties in computing the convolution of the fractional Green function and the forcing function [31,32].
Calculations of the transient vibrations of the analyzed beam with the fractional Kelvin-Voigt material are performed using the equation presented in the author's previous paper [32].
where l ck is the time constant.
This equation can be solved using the separation of variables method. Employing this method, the following equations of generalized coordinates are obtained D ð2Þ n n t ð Þ ð Þþx 2 n n n t ð Þ þ l ck D ðcÞ n n t ð Þ ð Þ The calculation of the general coordinates (Eqs (38)) for the beam using the fractional Kelvin-Voigt material is performed numerically using the same method as for the beam with the fractional Zener material. Thus the Eq. (38) is converted into a system of differential equations of a mixed integer and fractional order The solution of the system of equations above can be simplified using the identity Eq. (36) The systems of mixed integer and fractional order differential equations Eqs. (35) and (39) can be partitioned into two separated systems of equations solved simultaneously. The first system is of integer order differential equations while the second one is of a fractional order. Two different algorithms can be employed to solve each system.
The systems of Eqs. (35) and (39) are solved using a procedure for solving differential equations of a mixed integer and fractional order. The procedure is written by the author and implemented in the ''Matlab'' package. The fractional order differential equations are solved using the trapezoidal rule for the fractional Caputo derivative worked out by Diethelm et al [54]. The integer order equations are solved using the wellknown Adams-Bashforth-Moulton predictor-corrector method [55,56].
The roots of the characteristic beam equation (21) are calculated using the ''FindRoot'' procedure of the ''Mathematica'' package whereas the roots of the characteristic equation (24) are computed using the ''NSolve'' procedure of this package.

Investigation of the effect of the fractional derivative order on transient vibrations of a beam
The equations achieved in the previous section are used to study the impact of the order of the fractional derivative c on the dynamic properties of the analyzed beam. Studies are made of the Plexiglas beam with (20)). First, the impact of the order c of the fractional derivative on the roots of the characteristic equation was analyzed. These roots, namely s n ¼ Àr n þ iX n , were calculated for a ¼ 1:0, b ¼ 0:1, l c = 0.01841 s c and for two values of the coefficient b i.e. 0.0161 and 0.0175 s c . The calculations are performed for the first two vibration modes, which correspond to the first two roots of the Eq. (24). As mentioned before, these roots affect the damping coefficients and vibration amplitudes of a beam. For each mode, two complex conjugate roots are calculated.
The dependence of the calculated roots of the order of the fractional derivative c is shown in Figs. 2 and 3. The changes in the absolute value of the real part of the root are shown in Fig 2, where it can be seen that the absolute value of the real part of root, i.e. the damping coefficient, increases with an increasing order of the fractional derivative (Fig. 2). The changes in the value of the imaginary part of the root, which corresponds to the natural damped frequency X n , are shown in Fig 3. Note that the value of the damped natural frequency changes rather insignificantly with an increase in the order of the fractional derivative (Fig. 3).The damping properties of elements are often characterized by a dimensionless damping coefficient. The dimensionless damping coefficient f can be defined as the ratio of the real part (damping coefficient) to the imaginary part (damped natural frequency) of the root of the Eq. (24). The dependence of the dimensionless damping coefficient on the order of the fractional derivative for the first and second mode of vibrations is shown in Fig. 4. The dependence of the dimensionless damped natural frequency (X n =x 1 ) on the order of the fractional derivative c for two first vibration modes is presented in Fig. 5. It can be seen that the value of the first damped dimensionless natural frequency virtually does not dependent on the order of the fractional derivative.
In the next step, the effects of the order of the fractional derivative c on the transient responses of the beam are analyzed. The responses of the beam are calculated using a system of equations (Eq. (35)), and it is assumed that the distributed load qðx; tÞ ¼ 0. The beam is excited by its base motion. The motion of the beam base is determined by the following function of time where w 0 is the excitation amplitude, e is the angular acceleration.
The dimensionless beam deflection w=w 0 versus the dimensionless time parameter s is computed for the point at the end of the beam on which the mass element is located (Fig. 1)   As can be seen in Tables 1 and 2 , the maximum vibration amplitudes of the beam responses computed for b ¼ 0:0175 s c are greater then those computed for b ¼ 0:0161 s c . This is because for a constant value of the coefficient l c , the damping level depends on the value of coefficient b (see e.g. [44]), namely if b increases, the damping level decreases. Moreover, it is visible that the maximum amplitudes decrease with an increasing order of the fractional derivative. The results obtained show that with an increase in the order of the fractional derivative, the vibration amplitudes decrease (Tables 1 and 2, Figs. 6, 7, 8, and 9). This effect is more noticeable in the area of the second resonance, especially for lower values of the order of the fractional derivative. Moreover, it can be seen from the data in Tables 1 and 2     the amplitudes obtained for the material model formulated using integer order derivatives are significantly lower than those obtained using fractional derivatives.
In order to evaluate the differences between the calculations results obtained for the beams with the fractional Zener and Kelvin-Voigt models, additional calculations for a similar beam with the fractional Kelvin-Voigt model, described by Eq. (41), are performed: where l ck is the time constant. In these calculations, it is assumed that the loss tangent for the fractional Zener and Kelvin models is equal for the first undamped natural frequency of the beam under study. The loss tangent for the fractional Zener material is expressed as [44] tan dðmÞ where m is the frequency of consideration. The loss tangent of the fractional Kelvin-Voigt material for the frequency m can be calculated using equation [44] Therefore, the parameter l ck may be obtained from the expression below with m ¼ x 1 , where x 1 is the first undamped natural frequency of the beam, which is calculated using Eq. (22).
The time histories for the beam with the fractional Kelvin-Voigt model are calculated using Eqs. As can be seen from the Tables 1, 2, 1, and 4, the maximum vibration amplitudes of the first mode of vibration have virtually the same value for the results obtained with the fractional Zener and Kelvin-Voigt material models, which can be expected assuming that the loss tangent is equal for both materials for the first undamped natural frequency. It can also be seen that for the second mode of vibration, the maximum vibration amplitudes calculated for the fractional Kelvin-Voigt model are smaller than those for the fractional Zener model. As the order of the fractional derivative increases, the difference between the maximum amplitudes calculated using these models increases. These differences are greater for the coefficient b ¼ 0:0161 s c than for b ¼ 0:0175 s c , i.e. for greater material damping. In the case where b ¼ 0:0161 s c , the relative difference between the maximum amplitudes obtained using both models for c ¼ 0:4; 0:6 and b ¼ 0:05 is less than 3:1%, whereas for b ¼ 0:1 it is less than 9:7%. In the case where b ¼ 0:0175 s c , the relative difference between the maximum amplitudes calculated using aforementioned material models for c ¼ 0:4; 0:6 and b ¼ 0:05 is less than 0:9%, whereas for b ¼ 0:1 it is less than 0:8%.
From the time histories presented in Figs. 6, 7, 8, 9, 10, 11, 12, and 13 one can see that the corresponding time histories obtained for both material models are similar, and virtually identical for some physical parameters of the beam.
Thus, it can be stated that for some physical parameters, the use of of both material models generates very similar results. Therefore, the simpler fractional Kelvin-Voigt material model may be used instead of the fractional Zener model.  The calculations performed reveal that use of the fractional damping impacts on the time histories of the system. The calculated beam time histories show that for some values of the order of the fractional derivative c, the amplitudes occurring in the area of the second resonance are greater than those obtained in the area of the first resonance. In contrast, in the material models described by integer order derivatives, these amplitudes are considerably lower than those calculated in the area of the first resonance. A similar behavior of transient beam responses was demonstrated for the fractional Kelvin-Voigt material of the beam in the author's previous paper [32]. The occurrence of higher vibration amplitudes in the area of the second resonance and the achieved responses of the beam should be verified by experimental researches.
Additionally, studies are carried out to evaluate the difference between the results obtained for the calculations made using the fractional Zener model and the fractional Kelvin model. It is shown that for some physical beam parameters, the calculation results obtained using both models are virtually the same, meaning that the simpler, fractional Kelvin-Voigt material can be used instead of the fractional Zener material model.
Additionally, it has been found that the numerical computation time of the transient beam vibrations for the fractional Zener material model is much longer than that for similar dynamics computations using the Kelvin-Voigt model. This is because the order of the fractional differential equations is higher for the fractional Zener material model than for the fractional Kelvin-Voigt model. As a result, more fractionalorder differential equations are used for the computational method for the fractional Zener model than for the fractional Kevin model (see Eqs. (35) and (39)), which considerably increases the numerical computation time and complicates the algorithm for solving differential equations. In the author's opinion, this is a   Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

Appendix A
As is mentioned above, the equations of motion of the beam are derived by employing the Euler-Bernoulli beam theory , according to which rotation of element and shear deformation are disregarded. A free-body diagram corresponding to small element of length dx is shown in Fig. 14.
The force equation of motion in the vertical direction (axis Oz) gives Aq o 2 wðx; tÞ Canceling the appropriate terms, the equation above can be written as The moment equation of motion about the axis perpendicular to the plane Oxz, and passing through the centroid of the element dx (see Fig. 14) has the form Canceling the appropriate terms and disregarding terms having second powers in dx, the following relationship is derived Differentiating equation (48) Substituting for the shear force derivative (Eq. (46)) into equation (49), the following expression is derived o 2 Mðx; tÞ ox 2 ¼ Aq o 2 wðx; tÞ ot 2 À q x; t ð Þ ð50Þ Next, the relationship between the bending moment M(x, t) and the deflection w(x, t) should be obtained.
For an Euler-Bernoulli beam, the strain is expressed as Substituting relationship (51) into equation (1) and multiplying both sides by z and integrating over the cross-section area Next, substituting equation (50) into equation (54), the following relationship is obtained Appendix B A free-body diagram of a mass element attached to the free end of a beam is shown in Fig. 15. It is assumed that the mass center of the element coincides with the end of the beam, i.e. the distance BB 0 is neglected The rate of change of the linear momentum of the mass element m p in the vertical direction (axis Oz) gives m p D ð2Þ w B ðtÞ ð Þ¼ÀVðl; tÞ ð57Þ The rate of change of the angular momentum of the mass element m p with respect to its mass center B is where /ðx; tÞ is the rotation angle of the beam cross section, and for an Euler-Bernoulli beam equals the slope angle (see Fig. 15 The transverse translation of the element mass center w B ðtÞ is the sum of the rigid body and the relative translation [47], namely w B t ð Þ ¼ w l; t ð Þ ¼ w bs t ð Þ þ w rl l; t ð Þ ð62Þ Then, second derivative of w B t ð Þ with respect to time is