An Analysis of the Dynamical Behaviour of Systems with Fractional Damping for Mechanical Engineering Applications

: Fractional derivative models are widely used to easily characterise more complex damping behaviour than the viscous one, although the underlying properties are not trivial. Several studies about the mathematical properties can be found, but are usually far from the most daily applications. Thus, this paper studies the properties of structural systems whose damping is represented by a fractional model from the point of view of a mechanical engineer. First, a single-degree-of-freedom system with fractional damping is analysed. Speciﬁcally, the distribution of the poles and the dynamic response to several excitations is studied for different model parameter values highlighting dissimilarities from systems with conventional viscous damping. In fact, thanks to fractional models, particular behaviours are observed that cannot be reproduced by classical ones. Finally, the dynamics of a machine shaft supported by two bearings presenting fractional damping is analysed. The study is carried out by the Finite Element method, deriving in a system with symmetric matrices. Eigenvalues and eigenvectors are obtained by means of an iterative method, and the effect of damping is visualised on the mode shapes. In addition, the response to a perturbation is computed, revealing the inﬂuence of the model parameters on the resulting vibration.


Introduction
In the context of dynamics, fractional models allow us to easily represent the vibratory behaviour of elements that would otherwise require complex formulations, such as multielement or hereditary models, because they are able to reproduce correctly the damping mechanisms that come into play using a low number of variables [1,2]. Fractional models are specially advantageous for polymeric materials that show some level of dependence to frequency and arise naturally, for example, from certain motions of Newtonian fluids [3] or the molecular theories that predict the behaviour of certain type of polymeric materials [4]. In fact, fractional models are used to capture with more ease the viscoelastic nature of materials such as rubber or concrete [5] whose behaviour was previously modelled with a power law [6,7]; fractional operators appear in the non-linear stress-strain relation of metals [8]; and viscoelastic constitutive models based on fractional derivatives have been proposed to reproduce the time dependent behaviour of real materials [9][10][11][12][13][14].
Several works treating the behaviour of dynamic systems following fractional models from a mathematical perspective can be found in the literature. For instance, in [15], a fractionally damped single-degree-of-freedom (1 DOF) oscillator was analysed using the Laplace transform. It was found that a system of this kind exhibits nine distinct behaviour cases in opposition to the three shown by a traditional oscillator. Also, in the [16][17][18] series, the behaviour of a fractional oscillator whose derivative order was between one and two was studied, first when subjected to initial conditions and, after, when applying a sinusoidal force to the system. The main conclusion of the series was that a fractional oscillator like the one analysed presented damping even if not explicitly stated.
Computing methods have been widely proposed as well, mostly numerical [19] due to the difficulty to obtain closed-form solutions. Apart from the Laplace transform, which has been already mentioned, the Fourier transform [4], an eigenvector expansion in the state space [20], modal synthesis [21] or finite elements [22,23] have been used to obtain the response of systems presenting fractional damping.
In general, these types of systems are not analysed from the engineering point of view, even if fractional models have been used in the actual engineering practice, for instance, to linearise equations or to reduce the number of parameters needed to represent elastoplastic [24] or viscoelastic behaviours [11]. For this reason, the objective of this work is to shed light on the dynamic behaviour of a system whose damping is modelled using a fractional model. To this aim, first the theoretical background is presented with the aid of a 1 DOF system, for which the distribution of poles is studied and whose dynamic response to different excitations is computed. Instead of focusing on obtaining the analytical response for certain values, the behaviour of the system is analysed numerically in a range of values in order to identify trends. The goal is to identify the differences between a system whose damping is modelled as fractional and one following a traditional viscous damping formulation. This understanding is used afterwards to represent the dynamics of a machine shaft supported by two bearings with fractional damping.

Theoretical Background: A 1 DOF System with Fractional Damping
A single-degree-of-freedom mechanical system that presents fractional damping follows the equation of motion where m, c and k are respectively the mass, damping and stiffness of the system, F(t) is an external force and α is the order of the fractional derivative that, in order to describe a viscoelastic damper, is considered to be between zero and one-if the order of the derivative were zero, the damper would behave as a spring and would be considered "elastic"; if one, the damping would be viscous. The fractional derivative of a function f (t) is computed following the Caputo definition using the gamma function Γ as where n − 1 ≤ α < n. In this work, the Caputo definition is preferred to the Riemman-Liouville because the constant terms produced by the Laplace transform of the former have direct physical meaning [15] whereas the interpretation of such terms in the latter is more complex [12]. In the case under study, as the order of the derivative α is between zero and one and, in consequence, n = 1, the fractional derivative in (2) becomes For more detail about the definition or computation of fractional derivatives, any of the classical texts on the subject [2,[25][26][27] is recommended.

Poles of the System
One of the peculiarities of a system with fractional damping is the nature of its poles. In order to compute them, (1) is expressed in the Laplace domain, where adopts the form if F(t) = 0 is considered, X(s) = L {x(t)} being the Laplace transform of the displacement x(t).
Then, the characteristic equation is solved. Supposing that α is a rational number and can be therefore expressed as the quotient of two integers α = p/q, (5) becomes the polynomial of order 2q from which the poles of the system can be obtained. Resorting to the state space is a more elegant way to determine the poles of the system [8,20,28,29], but the presented process allows us to focus on the qualitative properties of the fractionally damped oscillator keeping the mathematical complexity to a minimum. From (6), one could think that a fractionally damped system has 2q poles but, in reality, only two of them-a complex conjugate pair [15], actually-are solutions to the original system (5), the rest are extraneous solutions due to the solving process.
The fact that (5) has only two complex conjugate poles can be explained from the physics of the real system it models: a 1 DOF system can only vibrate at a single frequency. Also, the fact that the poles of the system are complex conjugates implies that the response is always oscillatory and that there are no overdamped or critically damped cases as defined for a traditional 1 DOF system [20]. From a mathematical perspective, this phenomenon is related to the fact that the solutions of fractional order differential equations are Mittag-Leffler functions, whose properties explain the oscillatory dynamics of fractionally damped systems [16,29]. This feature will be further discussed when studying its dynamic response in Section 2.2.
In order to better understand the behaviour of the system, the variation in the position of the poles with the order of the derivative α is studied. Without loss of generality, for the numerical computations the mass m, damping c and stiffness k parameters are considered unitary (in any coherent system of units) while the value of the order of the derivative α ranges between zero and one, the completely elastic and the completely viscoelastic behaviours respectively.
If the poles of the system for different values of α are plotted in the complex plane, the shape shown in Figure 1 is obtained. In it, the solutions obtained from the variable change in (6) are shown in black while the ones that satisfy (5) are highlighted in red.
For every value of α, the centroid of each of the two distributions shown in Figure 1 computed as s j being the jth pole of the system and q the amount of solutions of (6) for each distribution, the extraneous ones included, fall in ±ω 0 = ± √ k/m. Another conclusion that can be drawn from the variation of the position of the poles with the value α is that, contrary to what could be expected, the number of decimal positions used to express α and, thus, the p and q values used to approximate it as a rational number, do not affect the position of the poles. This means that considering α or a slightly perturbed value α + δ provides the same results both in terms of poles and, in consequence, of time response. Then, the relationship of the damped frequency of the system, represented by the imaginary part of the poles, with the order of the derivative is studied. As shown in Figure 2, the damped frequency of the system ranges from when α = 0 (it should be noted that in this case, as the damper behaves like a spring, the system will be undamped) to the well known formula for the viscous case when α = 1 where ξ = c/(2 √ km) and ω 0 = √ k/m. The damped frequency decays with the order of the derivative α because, due to the viscoelastic nature of the model, the lower this order, a bigger part of the damping behaves as stiffness [30].
Then, the value of the damping parameter c is changed keeping the order of the fractional derivative constant, in order to study the evolution of the damped frequency ( Figure 3). In contrast to what happens in systems with viscous damping, increasing the damping parameter c in a fractionally damped system increases the vibration frequency. This can be explained again by the viscoelastic nature of fractional damping: while in a system with viscous damping a change in the value of the damping parameter c does not affect the stiffness; if the damping is fractional, both magnitudes are interrelated and an increase in the damping parameter c also stiffens the system. It is worth noting that, even if the tendency is the same for different values of the order of the derivative α-increasing the damping parameter c produces an increase in the damped frequency-the relationship between the two magnitudes is not linear and varies depending on the value of the order of the derivative.

Dynamic Response
Another aspect in which a system with fractional damping differs from one showing viscous damping is the dynamic response. These differences are analysed in the 1 DOF system, first, when it is subjected to initial conditions and, after, when a impulse or step force is applied.
In order to make the analysis extensible to systems with several degrees of freedom or subjected to any type of excitation for which analytical solutions are not available, the process to solve the differential equation (1) numerically is presented. The analytical response of a 1 DOF system to the cited excitations-initial conditions or impulse and step forces-can be found, for instance, in [2,20,29,31].
That said, the Grünwald-Letnikov method is used to compute the values of the fractional derivative for a given time t. Its objective is to take advantage of the definition of the derivative and approximate the gamma functions with some recursively computed coefficients. The derivative of order α is, thus, expressed as where S = t/∆t is the number of samples. Considering the A j coefficients can be computed recursively by 13) and the derivative of order α takes the simplified form which can be easily introduced in a numerical method, a forward Euler in this case.
It should be noted that fractional operators are convolution integrals and, as such, they have memory, making it is necessary to store the previous history so as to be able to compute the response in a certain step [1]. If needed, techniques as the one described in [30] could be used to reduce the storing requirements.
The presented methodology is followed in the next sections to compute the response of the fractionally damped 1 DOF system under study to initial conditions and to a impulse and step force, two of the most typical excitations.

Response to Initial Conditions
The response of a fractionally damping system whose m, c and k parameters are unitary and subjected to the initial conditions {x 0 , v 0 } = {1, 0} is shown in Figure 4. As predicted from the poles of the system, the response is oscillatory and decays in time. The main difference that a system with fractional damping presents in comparison to one that shows viscous damping regarding its response to initial conditions is the appearance of a non oscillatory term that decays in time [19]. This effect can be seen in Figure 5, where the response in the frequency domain, obtained by the Fast Fourier Transform of the time response, is represented for different values of the order of derivative α. During the initial part of the response (Figure 5a), only a peak at the vibration frequency can be noticed; in the final part of the decay (Figure 5b), however, it is possible to identify also a peak in the spectrum at zero frequency, that is related to the non oscillatory term, together with the peak at the vibration frequency. From a physical point of view, the appearance of a non oscillatory term has two implications. First, it means that the equilibrium position of the system is not constant but it varies with time. Taking into account that the equilibrium position is defined by the mass and the stiffness of the system and, in this case, the mass is constant, the single explanation to this phenomenon is that the stiffness of the system varies with time. Secondly, it implies that the system does not oscillate with decreasing amplitude until it finally stops, but that it reaches a point in which it approaches the original equilibrium position without oscillating.
Another aspect in which a system with fractional damping differs from one with viscous damping is in the fact that a critical damping ratio in which the system does not oscillate around the equilibrium position cannot be defined [20,30]. The response always crosses the origin of coordinates at least once due to the nature of the poles of the system, that are always complex conjugates. For this reason, fractional damping allows us to model dampers whose response only crosses the equilibrium position a single time, a case that cannot be represented with accuracy using a traditional model. For example, the 1 DOF system under study presents this behaviour when the damping parameter c = 2 and the order of the derivative α = 0.9 ( Figure 6).

Response to Impulse
The response to an impulse is of interest because it provides a way to define critical damping. According to [20], critical damping can be defined for systems with fractional damping as the value for which the response is tangent to the axis of zero displacement when an impulse force is applied to the system. This means that for values of damping greater than this the response does not change sign, even if it still oscillates, around a position that varies with time in this case (Figure 7). Following this reasoning, the value of critical damping can be estimated for different values of the order of the derivative α. For example, in Figure 8, the response of the 1 DOF system with unitary mass and stiffness parameters and critical damping when subjected to a impulse is presented. For low values of α the response always changes sign no matter how high the added damping: increasing the damping parameter makes the system stiffer and the effect of the increase in the natural frequency is greater than that of the dissipation.

Response to Step
The response of the system with fractional damping to a step force has the peculiarity that it stays under the response of a system with viscous damping. As an example, the response of the 1 DOF of freedom system with unitary mass, damping and stiffness constants under study is presented in Figure 9 depending of the order of the derivative α for a unit step excitation. It can be noted that for systems with fractional damping the overshoot decreases but the settling time is much longer, specially when decreasing the value of the order of the derivative.

Application: Bearing Support
In order to illustrate the properties of the fractional models presented in the previous sections in a real life application, a shaft supported by one bearing in each end is considered. The shaft is modelled by means of beam finite elements; both supports, following the typical model of a bearing [32], as a combination of a spring and fractional damper ( Figure 10). For the numerical application, a shaft of length L = 1 m and circular section of diameter d = 0.4 m made of steel with a Young's modulus of E = 210 GPa and density ρ = 7900 kg/m 3 is considered. The natural frequencies, mode shapes and response of the shaft when changing the parameters of the bearings are studied, the two bearings being identical (k 1 = k N−1 = k, c 1 = c N−1 = c and α 1 = α N−1 = α).

Natural Frequencies and Mode Shapes
The eigenpairs of a the shaft supported by bearings that present fractional damping satisfy the relationship where λ * r and φ * r are respectively the eigenvalue and the eigenvector of the rth mode, M the mass matrix of the system, that corresponds to the mass matrix of the shaft as the supports are considered massless, K shaft and K * supp (ω) the stiffness matrices of the shaft and the bearing supports. This last complex matrix K * supp (ω) is frequency dependent and is defined as where B holds the positions of the bearing supports, in this case the first and penultimate (N-1) degrees of freedoms, the ones that represent the vertical displacement of the first and last nodes, respectively and e j stands for the unit vector that corresponds to the jth degree of freedom. As the matrix K supp is frequency dependent, an iterative method is needed to obtain the eigenpairs of the system: the reader can find the details in [33,34].
With the aim of studying the evolution of the vibration frequencies of the system with the order of the derivative α, the values of the stiffness and damping parameters of the bearing supports are set to k = 10 4 N/m and c = 50 N s α /m respectively. When the shaft is discretised with 60 beam elements in its length the results shown in Figure 11 are obtained. It can be observed that, as deduced for the 1 DOF system, changing the order of the derivative affects the stiffness of the supports and, thus, the vibration frequencies. Another aspect that should be taken into account regarding the eigenpairs is that the modes are complex due to the non-proportionality of the damping. When the value of the damping parameter c is low this effect is barely noticeable so, for the sake of clarity, its value is increased to c = 1000 N s α /m in the computation of the mode shapes presented in Figure 12. The phase of the mode shape in each point of the shaft is represented with the aim of bringing to light the complex nature of the modes: if the mode were normal, the phase would change abruptly from −π to π (or vice versa) in the nodal points, but it is not what occurs.
It should be noted also that the mode shapes are not affected by the order of the derivative α but the complex nature of the modes is more acute for high values of α as the effect of the damping is higher.

Dynamic Response
Next, the time response of the shaft to a perturbation from the base is computed following the procedure in [35] by means of a M-C-K Newmark method. The details of the implementation of such method can be found in [22], but, in short, it consists of expressing the equation of motion of the system in the instant n + 1 as and introducing the notation in order to use the Grünwald-Letnikov definition of the fractional derivative in the computations. In the case under study the mass matrix M is that of the shaft; the stiffness matrix K is the sum of the one of the shaft and the one of the supports where and the damping matrix C is given by the damping of the bearing supports as The time response of the supported shaft to a seismic excitation is computed as follows: first, the static response to a displacement of the rightmost end is computed; then, it is introduced into the presented Newmark scheme as an initial deformation u 0 . Four distinct cases are studied: The four cases are simulated for a time T = 1 s using 1000 samples and discretising the shaft with 60 beam elements. In order to analyse the difference between the cases, the displacement of the right end, to which a vertical displacement of 3 mm is imposed in the first instant, is plotted over time ( Figure 13). The first case represents the classical vibration damped in time: this behaviour can be modelled both with a viscous and a fractional damper. The motion of the second case, instead, is characteristic of a overdamped system and only achievable if the damping is viscous: a fractionally damped system will always oscillate, even if the response does not have to necessarily change sign, as the third case makes clear. Finally, the fourth case represents a more general situation in which the response consists of a combination of several modes and the rigid body motion. In this case, by tuning the stiffness and damping parameters and the order of the derivative, different behaviours can be reproduced.
The evolution of the shaft in time can be seen in the four videos provided as Supplementary Materials with this article.

Conclusions
This work discusses the behaviour of a 1 DOF system with fractional damping from the point of view of dynamics, specially focusing in the divergence of such a model of damping from a classical viscous one.
The first difference between a system with viscous and a fractional damping is that the poles of the latter are always complex conjugates, which implies that the response is always oscillatory and that, therefore, the concept of critical damping as defined for viscous systems is not applicable. Also, due to the viscoelastic nature of fractional damping, the damped frequency of the system decreases with the order of the derivative, as the system resembles more a pure viscous system, but increases with the damping parameter, as the contribution to the stiffness of the system increases too.
Several particularities can be observed regarding the time response as well. First, the response to initial conditions is characterised by the appearance of a non oscillatory term, which can be understood as a change in the position of the equilibrium position in time. Secondly, the response to an impulse allows us to define critical damping for this kind of systems as the amount of damping needed for the response not to change its sign. Finally, the response to a step function is characterised by its reduced overshoot but longer settling time is comparison to a system with viscous damping.
Despite their peculiarities, fractional models are useful to ease the modelling process and replicate behaviours that would otherwise require complex formulations. This is the case, for example, of the bearing support represented in this work with a spring and fractional damper for which behaviours that range from a traditional viscous case to an oscillating case with variable equilibrium position can be represented with a single model.
As a final remark, we would like to warn against using fractional models as mathematical artifacts that allow us to simplify the modelling process. It is advisable to check if the behaviour of the mechanical system under study follows the same trends as the assumed model, or at least, that the effects attributable to the model do not interfere in the computations.