Hysteretically Symmetrical Evolution of Elastomers-Based Vibration Isolators within α -Fractional Nonlinear Computational Dynamics

: This study deals with computational analysis of vibration isolators’ behavior, using the fractional-order di ﬀ erential equations (FDE). Numerical investigations regarding the inﬂuences of α -fractional derivatives have been mainly focused on the dissipative component within the di ﬀ erential constitutive equation of rheological model. Two classical models were considered, Voigt-Kelvin and Van der Pol, in order to develop analyses both on linear and nonlinear formulations. The aim of this research is to evaluate the operational capability, provided by the α -fractional derivatives within the viscous component of certain rheological model, to enable an accurate response regarding the realistic behavior of elastomeric-based vibration isolators. The hysteretic response followed, which has to be able to assure the symmetry of dynamic evolution under external loads, and at the same time, properly providing dissipative and conservative characteristics in respect of the results of experimental investigations. Computational analysis was performed for di ﬀ erent values of α -fractional order, also taking into account the integer value, in order to facilitate the comparison between the responses. The results have shown the serviceable capability of the α -fractional damping component to emulate, both a real dissipative behavior, and a virtual conservative characteristic, into a unitary way, only by tuning the α -order. At the same time, the fractional derivative models are able to preserve the symmetry of hysteretic behavior, comparatively, e.g., with rational-power nonlinear models. Thereby, the proposed models are accurately able to simulate speciﬁc behavioral aspects of rubber-like elastomers-based vibration isolators, to the experiments.


Introduction
Fractional calculus (FC) studies the differentiation and integration of an arbitrary real or complex order of the differential operator. The starting point of FC's use was a discussion between Leibniz and L'Hospital (1695) concerning the calculation and significance of the derivative of the power function, which also became the name of FC. FC has been largely applied to pure mathematics and its application in physics or engineering has remained extremely limited until, by association with the theory of chaos or fractals, it has increased the interest in this branch of analysis. The reputation enjoyed by FC is still exotic, and for many researchers, it is still unclear what its practical use is. Fractional derivatives take into account the previous evolution of the variable, unlike the standard differential operator, which is limited to taking into account local time history. This feature allows for a better description of system dynamics [1].
Taking the previous aspects into account, the next paragraphs contain a briefly survey of some relevant works within the areas of both computational aspects, and some practical applications on vibration isolation analysis, of FC. Thus, in article [1], there is presented a dynamic statistical analysis of a system composed of a large number of masses. The conclusion of this article shows that, while the individual dynamics of each material point has an integer-order behavior, global dynamics reveal both the existence of the integer dynamics and the fractional dynamics. Babiarz and Legowski, in the paper [2], developed a fractional model of the human arm, and is an approach representing an attempt to consider the properties of damping in the simplest form.
Inside the editorial paper [3], Zhou et al. had proposed a selection of papers presenting many applications of fractional dynamics in different research domains. Three areas are distinguishable: Dynamical analysis (leading to fractional differential equations and systems), fractional-order control theory, and applications considering models with fractional calculus.
The work [4] presents some ways to use Matlab software with the associated numerical methods to simulate and compute the fractional derivatives. These methods are based on series technique (Fourier and Taylor in our case). In addition, some numerical examples were presented. The fractional calculus offers a number of axioms and methods to extend the derivative notions from integer n to arbitrary order α, under favorable conditions [5]. In paper [5], the author has intended to solve a linear homogeneous differential equation with constant coefficients by means of N-fractional calculus method, and the results show that singular differential equations can be solved by means of the N-fractional calculus method.
Koh and Kelly, in Reference [6], are studying a mechanical oscillator consisting of a concentrated mass and a Kelvin-type element with fractional behavior, modeling an elastomeric bearing in groundwork isolation systems. For a dynamic analysis of a fractional oscillator developed efficient numerical multi-step schemes in the time domain, in good agreement with the Laplace and Fourier solutions. The model with fractional derivative agrees, to a satisfactory extent, with the experimental results.
Freundlich, in paper [7], is studying the vibration of a simply supported beam with a viscoelastic behavior of fractional order. The beam, considering the Bernoulli-Euler hypotheses, was excited by the supports movement. In this model, the Riemann-Liouville fractional derivative (RLFD) of order 0 < α ≤ 1 was used. Firstly, the steady-state vibrations were analyzed and therefore the RLFD with lower terminal at -∞ was considered, in order to simplify solution of the FDE and to enable direct computation of the amplitude-frequency characteristic. Final results showed that the appropriate selection of both damping coefficients, and respectively, of fractional derivative order, offers more accurate fitting of the dynamic characteristic comparatively with the model using the integer order derivative.
The fractional equations of the mass-spring-damper system were considered by Gómez-Aguilar et al. in Reference [8]. The results had shown viscoelastic behaviors of the mechanical components (producing, at different scales, temporal fractality). This proves the material heterogeneities in the mechanical components.
Use of the fractional order models in bioengineering applications is presented in paper [9], where a method to obtain a numerical solution, considering the power series expansion was provided. Fractional dynamic of such system leads to different equations and the applications are solved with Matlab/Simulink software. The study offers several illustrative cases that can be used not only in bioengineering, but so too in other disciplines, where it is possible to apply the fractional calculus. In addition, with the previous mentioned paper [9], within chapter [10], Petráš had developed methods to obtain the fractional derivative, fractional integral solution of a fractional differential equations using Matlab. Illustrative examples, together the correspondent Matlab code sources, were also presented.
The researches of Rossikhin et al. dealt with dynamic behavior of mechanical oscillators (linear and nonlinear)-in [11], respectively, with force-driven vibrations of nonlinear oscillators-in Reference [12], where the constitutive equations use fractional derivatives. The authors had concluded that using fractional derivatives is possible to do an approximate analysis of the vibratory regimes of the oscillator.
The nonlinear response of a plate endowed loaded with a random force using fractional derivative is presented in Reference [13] by Malara and Spanos. A statistical linearization is used to determine an approximate solution of the governing equations of the plate vibration. Eigenvalues analysis offers a time-dependent representation of the response. In this way, it is possible to obtain the set of nonlinear fractional differential equation. The linearization is considered in a mean square sense.
Inside study [14], the stochastic response of a structural systems has been studied, with a single-degree-of-freedom. The excitation was stationary and non-stationary. The damping is considered with fractional behavior. The final equation of motion becomes a set of coupled differential linear equations, involving more degrees of freedom. The number of the additional degrees of freedom depends on the discretization of the fractional derivative operator. Using these methods, the stochastic analysis becomes straightforward and simple. It is mentioned that, for the most current engineering interest, it has been proven that the second-order statistics response can be obtained in a closed form, in a simple manner.
In paper [15], Cajić et al. makes an analysis of a nanobeam free damped transverse vibration. The method used is the non-local theory of Eringen and fractional derivative model for viscoelasticity. The classic Euler-Bernoulli beam theory and constitutive equation using fractional order derivatives that offer the form of the motion equation of a nanobeam. Different values of fractional order parameter are used in order to obtain the time dependent behavior.
The actual literature also provides a suite of theoretical approaches related to fractional differential equations. Hereby, Ma [16] proposed the Adomian decomposition method and a computational technique for the study of the impulsive fractional differential equations. Continuing this research in Reference [17], to obtain an exact solution for nonlinear fractional problems, the authors combined the double Laplace transforms with the Adomian decomposition method. The new method is applied to solve regular and singular conformable fractional coupled Burgers' equations, and illustrated the effectiveness of the proposed method with some examples.
Finally, in this section, and in addition to previously mentioned works, comprehensive information regarding both theoretical and computational aspects of FC can be found within the texts [18][19][20][21][22]. These works also contain some applications of FC on various important areas of engineering and physics.
Linear and nonlinear models regarding the elastomers-based vibration insulation devices can be found in papers [23][24][25], where they had earnestly presented the influence of dissipative characteristic on isolator dynamic behavior. In addition, the research trials, as well as the efforts of the author, within the area of vibration isolation using the passive systems, was presented by the papers [26][27][28][29][30], which contain both computational and experimental results, according with proposed linear and nonlinear approaches.
Usually, numerical investigations have a sensible aspect related to the required computational resources, optimal balancing between hardware performances, and computing time. Most of the researchers made additional searches into computing time optimization, beyond that of the mainly numerical procedures development [31]. Regarding this aspect and within the present study, the authors noted that the computational cost required by evaluation of α-fractional derivative increase the total simulation time to a certain measure, directly depending on both FC formulation and solving the numerical approach had adopted within the computational model. Computational cost optimization was not framed by the aim of this study, but nevertheless, some related-to investigations done by the authors show a relatively small increasing of total simulation time for FC-based algorithm comparatively with the linear case. Taking into account an average 6000 points single dimension grid for each simulation, it has to be noted that the FC had required 1.12...1.23 s, in the same time with 0.88...0.93 s by the classical algorithm (hardware and numerical FC formulation were kept identical).
This paper was structured based on the next sections as follows: Theoretical basics of α-fractional derivative calculus and classical models related on FDE's within vibratory systems analysis, computational assessments regarding fractional-order dissipation, discussions, and final concluding remarks.

Materials and Methods
The basic approach of fractional calculus consists by the Riemann-Liouville fractional integral (RLFI) [4,9,10] with α ∈ (0, 1] and t > a, where a is the initial value (usually a = 0). Associated derivative is able to result through the Lagrange rule for differential operators. Evaluating the nth order derivative over the (n − α) order integral, yields the α order derivative (related to RLFD) as follows with α ∈ (n − 1, n], α > 0, notation D denoting d/dt and Γ(•) means the Gamma function, defined as follows Taking into account the initial value a set to 0, that denotes the most common case, the RLFD, Equation (2), becomes [2] where in it was assumed the restriction of derivative order α ∈ (0, 1]. Following the RLFI previous expressions, a commonly compact definition, used for fractional derivatives, is [6] Caputo, in 1990's, supposed another way to define the fractional derivatives, actually known as Caputo Fractional Derivative (CFD) [4,8] where α = β + m, β ∈ (0, 1] and m denotes the nearest integer less than α. The Caputo definition does not require the fractional order initial conditions to define, and thus, it is often preferred in solving differential equations. A.K. Grünwald and A.V. Letnikov propose a new generalization of integer derivative by approaching the integer difference of the fractional case [4]. Supposing the binomial coefficients as an iterative formulation [10] x y = c (x) and considering the real valued binomial coefficients based on Gamma function [10,16] because that allows a generalization to non-integer arguments the Grünwald-Letnikov Fractional Derivative (GLFD) is [4,9,10] with the notation [•] meaning the integer part of the argument. Obviously, y(t) = 0 for t < 0 within transient problems, thus −∞ D α t y(t) becomes identical to 0 D α t y(t), [6]. For analyses within the area of vibration isolation, into the context of the viscoelasticity theory, the fractional order α, known as memory parameter, was assumed to have a value varying from zero to one.
The fractional model, used within computational analyses, followed the conventional Voigt-Kelvin rheological schematization. The constitutive equation, embedding fractional order formulation into dissipative term, were assumed to have the following expression [30] with damping constant c > 0, G o denoting the shear modulus, and τ, γ meaning the stress and the shear strain, respectively. It has to be noted that the authors of the paper [6], using the less-square fit procedure for a set of experimental data, showed that the fractional Voigt-Kelvin model fits the data much better than the hysteretic, classical Kelvin and standard linear solid models. Three practical applications of FDE for mechanical vibration analyses [30], embedding fractional order terms, are briefly presented as follows. The first case is the fractional order Van der Pol oscillator, having both inertial and dissipative α-fractional terms or supposing only the fractional damping ..
Second case relates to the equations of fractional order Duffing-like systems, with a dissipative term or without damping The last exemplification case presents an interesting system with small fractional and delayed damping, modeled by the following expression ..
In Equations (11)- (15), having mass normalized formulations, δ denotes damping coefficient, a > 0 is the control parameter of nonlinearity, b relates to squared pulsation of equivalent linear system, and f (t) is the excitation term. The parameters β i are the control coefficients related to each damping component. It has to be noted that extensive examples of FDE can be found in the works [11,12] and within theirs bibliography. One challenge, especially for practical simulations, consists of settling the correct dimensionality of the differential equations. In order to supply this, an additional parameter µ that must be included within the fractional temporal operator was proposed in the work [8] where parameter µ has a dimension of seconds. According to Reference [8], parameter µ relates to the temporal components within the system (those components implying changes to the time constant of system). Clearly, taking µ = 1 into Equation (16) yields ordinary integer temporal derivative operators.
According to Reference [8], µ can be estimated by following expressions for a mass-spring ensemble, for a damper-spring ensemble, and for a complete mass-spring-damper ensemble, with k, m denoting the rigidity and mass, respectively, and c meaning the damping coefficient.
In the view of this papers goal and taking into account the previous assertions regarding α-fractional derivatives, a α-fractional derivative SDoF system based on Voigt-Kelvin linear schematization, can be modeled with the expression [30] 1 where the parameters have aforementioned significations.
Following the aim of this study, to characterize the hysteretic behavior of vibration isolator based on α-fractional derivative approaches, the authors consider a differential formulation related to Equation (20), but ignore the inertial term. This hypothesis were justified by the assertions in paper [11], and was also mentioned and detailed by the analyses in paper [30], where an extended investigation regarding the implications of α-fractional derivative within the classic linear was developed, Duffing and Van der Pol oscillators, respectively. Different values for the main parameters were adopted and the responses, in terms of time, spectral, and hysteretic characteristics were discussed.
Into this paper, the symmetrical trends in evolution of a vibration isolator subjected to harmonic excitation were studied. The capability of this α-fractional derivative approach to simulate both dissipative and conservative character of an isolator through a singular term was investigated, involving only the damping, but α-fractional-typed. The analyses were focused on evaluation of the equivalence between α-fractional-derivative and pseudo-linear systems, but underlined the maintaining of symmetry of hysteretic behavior. This last aspect is very important because of the practical observations (including instrumental investigations), which highlight the symmetry in hysteresis loops independently by dynamic excitation type or material characteristics. In this view, the α-fractional approach comes with the advantage of a unitary expression that intrinsically embeds damping and stiffness parameters, but keeping the realistic behavior (in term of hysteresis symmetry).

Results
As was previously mentioned, a simple dissipative oscillator was adopted, according to classical Voigt-Kelvin schematization, Equation (20), with α-fractional derivative in damping term c · 0 D α t x(t) and null stiffness. The inertial term was ignored, the excitation was supposed to be kinematical, in terms of harmonic displacement x(t) with 1 Hz frequency, and the damping coefficient was unitary valued (c = 1). Computational investigations were conducted for order α ∈ {0.1, 0.3, 0.5, 0.7, 0.9}, and in addition, the ordinary case (α = 1) was considered. Fractional derivatives were evaluated using a GLFD-based algorithm.
Timed evolution and spectral composition of the response function were, respectively, presented in Figures 1 and 2. The results related to the regular derivative case (α = 1) were differently marked onto the graphs in order to facilitate comparative analysis. For spectral distribution, both magnitude and phase components because some changes in frequency evolution were considered and have more visibility into the phase diagram than the magnitude.

Results
As was previously mentioned, a simple dissipative oscillator was adopted, according to classical Voigt-Kelvin schematization, Equation (20), with α-fractional derivative in damping term  and null stiffness. The inertial term was ignored, the excitation was supposed to be kinematical, in terms of harmonic displacement x(t) with 1 Hz frequency, and the damping coefficient was unitary valued (c = 1). Computational investigations were conducted for order , and in addition, the ordinary case (α = 1) was considered. Fractional derivatives were evaluated using a GLFD-based algorithm.
Timed evolution and spectral composition of the response function were, respectively, presented in Figures 1 and 2. The results related to the regular derivative case (α = 1) were differently marked onto the graphs in order to facilitate comparative analysis. For spectral distribution, both magnitude and phase components because some changes in frequency evolution were considered and have more visibility into the phase diagram than the magnitude.   Figure 3 presents the response function in respect of the excitation displacement. The presence of the hysteretic loop can be observed, which acquires changes as a function of fractional order. In order to evaluate the linear-like evolution of the system, a Poincare map, with a snapshot increment identical to the time period of excitation signal, was overlapped on diagrams in Figure 3. Null scattering of dots related to Poincare maps highlight the correlated evolution of the system with the excitation frequency. It has to be mentioned that, for the investigations presented into this paper, the evaluation of Poincare maps is justified by the fact that the response function is identical with the system velocity, because the model only incorporates the unitary gained dissipative term (in fact, the response function in respect to the displacement is similar with the phase plane diagram).

Results
As was previously mentioned, a simple dissipative oscillator was adopted, according to classical Voigt-Kelvin schematization, Equation (20), with α-fractional derivative in damping term  and null stiffness. The inertial term was ignored, the excitation was supposed to be kinematical, in terms of harmonic displacement x(t) with 1 Hz frequency, and the damping coefficient was unitary valued (c = 1). Computational investigations were conducted for order , and in addition, the ordinary case (α = 1) was considered. Fractional derivatives were evaluated using a GLFD-based algorithm.
Timed evolution and spectral composition of the response function were, respectively, presented in Figures 1 and 2. The results related to the regular derivative case (α = 1) were differently marked onto the graphs in order to facilitate comparative analysis. For spectral distribution, both magnitude and phase components because some changes in frequency evolution were considered and have more visibility into the phase diagram than the magnitude.   Figure 3 presents the response function in respect of the excitation displacement. The presence of the hysteretic loop can be observed, which acquires changes as a function of fractional order. In order to evaluate the linear-like evolution of the system, a Poincare map, with a snapshot increment identical to the time period of excitation signal, was overlapped on diagrams in Figure 3. Null scattering of dots related to Poincare maps highlight the correlated evolution of the system with the excitation frequency. It has to be mentioned that, for the investigations presented into this paper, the evaluation of Poincare maps is justified by the fact that the response function is identical with the system velocity, because the model only incorporates the unitary gained dissipative term (in fact, the response function in respect to the displacement is similar with the phase plane diagram).  Figure 3 presents the response function in respect of the excitation displacement. The presence of the hysteretic loop can be observed, which acquires changes as a function of fractional order. In order to evaluate the linear-like evolution of the system, a Poincare map, with a snapshot increment identical to the time period of excitation signal, was overlapped on diagrams in Figure 3. Null scattering of dots related to Poincare maps highlight the correlated evolution of the system with the excitation frequency. It has to be mentioned that, for the investigations presented into this paper, the evaluation of Poincare maps is justified by the fact that the response function is identical with the system velocity, because the model only incorporates the unitary gained dissipative term (in fact, the response function in respect to the displacement is similar with the phase plane diagram).  The changes in slope of each hysteresis in Figure 3 indicate that, in the same time with increasing of α order, the system provide both dissipative and conservative characteristics, which quantitatively depend on the effective α value. This fact leads to the necessity of the pseudo-linear equivalent system evaluation, in terms of stiffness and damping of Voigt-Kelvin schematization. Taking into account the proportionality between the dissipated energy and the hysteresis loop area, it results in the equivalent-damping coefficient. In order to evaluate the stiffness, the hysteresis axis between the extreme horizontal coordinates was considered. For the loops in Figure 3, the correspondent axes were depicted in Figure 4. Analyzing the slopes of axes in Figure 4 clearly result in that α-fractional derivatives provide an additionally elastic component within system behavior. In the same time, the symmetry of the hysteresis loops through both the equality of the extremely values, and the common point of all axes on origin (0, 0) is evident.  The changes in slope of each hysteresis in Figure 3 indicate that, in the same time with increasing of α order, the system provide both dissipative and conservative characteristics, which quantitatively depend on the effective α value. This fact leads to the necessity of the pseudo-linear equivalent system evaluation, in terms of stiffness and damping of Voigt-Kelvin schematization. Taking into account the proportionality between the dissipated energy and the hysteresis loop area, it results in the equivalent-damping coefficient. In order to evaluate the stiffness, the hysteresis axis between the extreme horizontal coordinates was considered. For the loops in Figure 3, the correspondent axes were depicted in Figure 4. Analyzing the slopes of axes in Figure 4 clearly result in that α-fractional derivatives provide an additionally elastic component within system behavior. In the same time, the symmetry of the hysteresis loops through both the equality of the extremely values, and the common point of all axes on origin (0, 0) is evident.  The changes in slope of each hysteresis in Figure 3 indicate that, in the same time with increasing of α order, the system provide both dissipative and conservative characteristics, which quantitatively depend on the effective α value. This fact leads to the necessity of the pseudo-linear equivalent system evaluation, in terms of stiffness and damping of Voigt-Kelvin schematization. Taking into account the proportionality between the dissipated energy and the hysteresis loop area, it results in the equivalent-damping coefficient. In order to evaluate the stiffness, the hysteresis axis between the extreme horizontal coordinates was considered. For the loops in Figure 3, the correspondent axes were depicted in Figure 4. Analyzing the slopes of axes in Figure 4 clearly result in that α-fractional derivatives provide an additionally elastic component within system behavior. In the same time, the symmetry of the hysteresis loops through both the equality of the extremely values, and the common point of all axes on origin (0, 0) is evident.  At this point, it is evident that the α-fractional-based dissipative model is able to provide a behavior equivalent to a linear Voigt-Kelvin model, and vice-versa (see Figure 6). In fact, the fractional dissipative model is able to emulate a linear visco-elastic model, with direct linkage between α order and the pair of viscous and elastic parameters, respectively.

Diagram in
c 0 Figure 6. Schematic diagram related to the equivalence between α-fractional dissipative model and emulated linear Voigt-Kelvin model.
A computational application was developed for identification and evaluation of the proper values of the pseudo-linear system parameters. Different values for α-order between zero and one were considered, and both fractional and equivalent systems responses have been compared. The results indicate that the estimated relative error tends to zero, obviously depending by the computation precision. The diagrams in Figure 7 present two cases related to α-order valued to 0.4 and 0.7, respectively. It can be observed that both plots perfectly overlapped. At this point, it is evident that the α-fractional-based dissipative model is able to provide a behavior equivalent to a linear Voigt-Kelvin model, and vice-versa (see Figure 6). In fact, the fractional dissipative model is able to emulate a linear visco-elastic model, with direct linkage between α order and the pair of viscous and elastic parameters, respectively. At this point, it is evident that the α-fractional-based dissipative model is able to provide a behavior equivalent to a linear Voigt-Kelvin model, and vice-versa (see Figure 6). In fact, the fractional dissipative model is able to emulate a linear visco-elastic model, with direct linkage between α order and the pair of viscous and elastic parameters, respectively.
c 0 Figure 6. Schematic diagram related to the equivalence between α-fractional dissipative model and emulated linear Voigt-Kelvin model.
A computational application was developed for identification and evaluation of the proper values of the pseudo-linear system parameters. Different values for α-order between zero and one were considered, and both fractional and equivalent systems responses have been compared. The results indicate that the estimated relative error tends to zero, obviously depending by the computation precision. The diagrams in Figure 7 present two cases related to α-order valued to 0.4 and 0.7, respectively. It can be observed that both plots perfectly overlapped. A computational application was developed for identification and evaluation of the proper values of the pseudo-linear system parameters. Different values for α-order between zero and one were considered, and both fractional and equivalent systems responses have been compared. The results indicate that the estimated relative error tends to zero, obviously depending by the computation precision. The diagrams in Figure 7 present two cases related to α-order valued to 0.4 and 0.7, respectively. It can be observed that both plots perfectly overlapped. Taking into account aforementioned observation, it has to identify the correlation between α-order and the parameters of equivalent linear system. This analysis was performed by considering the equivalent stiffness and damping, respectively, and the results were depicted in Figure 8. In addition, the losses coefficient, as the ratio between dissipative and conservative parameters, was evaluated. The evolution of losses coefficient in respect to α-order was presented in Figure 9, where semi-log representation was adopted taking into account the relative large domain of function variability. It is evident that the losses coefficient g → +∞ for α = 1 and g = 0 for α = 0.   Taking into account aforementioned observation, it has to identify the correlation between α-order and the parameters of equivalent linear system. This analysis was performed by considering the equivalent stiffness and damping, respectively, and the results were depicted in Figure 8. In addition, the losses coefficient, as the ratio between dissipative and conservative parameters, was evaluated. The evolution of losses coefficient in respect to α-order was presented in Figure 9, where semi-log representation was adopted taking into account the relative large domain of function variability. It is evident that the losses coefficient g → +∞ for α = 1 and g = 0 for α = 0. Taking into account aforementioned observation, it has to identify the correlation between α-order and the parameters of equivalent linear system. This analysis was performed by considering the equivalent stiffness and damping, respectively, and the results were depicted in Figure 8. In addition, the losses coefficient, as the ratio between dissipative and conservative parameters, was evaluated. The evolution of losses coefficient in respect to α-order was presented in Figure 9, where semi-log representation was adopted taking into account the relative large domain of function variability. It is evident that the losses coefficient g → +∞ for α = 1 and g = 0 for α = 0.   Taking into account aforementioned observation, it has to identify the correlation between α-order and the parameters of equivalent linear system. This analysis was performed by considering the equivalent stiffness and damping, respectively, and the results were depicted in Figure 8. In addition, the losses coefficient, as the ratio between dissipative and conservative parameters, was evaluated. The evolution of losses coefficient in respect to α-order was presented in Figure 9, where semi-log representation was adopted taking into account the relative large domain of function variability. It is evident that the losses coefficient g → +∞ for α = 1 and g = 0 for α = 0.   An extended analysis was developed by taking into account the classical Van der Pol (VdP) expression for dissipative component, obviously supposing the α-fractional order derivative. The effects of fractional derivative in model behavior was followed, related to the observations and acquired for the previous simplest pseudo-linear model. This analysis was adopted based on the previous work evaluations [30], where it was considered for both Duffing and VdP non-linear models. The dissipative term within Duffing-type model enables only linear dependence in respect to the displacement derivative thus that the dynamic behavior cannot be differently influenced by the fractional derivative than the Voigt-Kelvin-based model. Hereby, in the frame of this study and taking into account the nonlinear damping expression, becomes interesting to analyze a VdP-type model, embedding fractional order derivative, related to its capability that emulate a classical VdP model. The following expression of VdP nonlinear model including α-fractional derivative was proposed for computational investigations [30] where f (t) denotes the response of the system, x(t) is the instantaneously displacement, and c 0 , k 0 denote damping and stiffness coefficients, respectively. In respect of the hypothesis within previous case analysis, the parameters c 0 , k 0 were identically valued as follows: c 0 = 1, k 0 = 0. Taking into account the nonlinear term in parenthesis within the left side hand of Equation (22), the analysis has to be conducted related to displacement magnitude. In this study, random values for magnitude were taken. The results presented correspond to the 3.3, 1.0, and 0.3 valued displacement magnitudes. The monitored parameters, for each case, were similarly adopted than the previous linear case: Timed evolution, spectral composition, hysteretic loops with Poincare maps overlapped, and hysteresis axis according to the pseudo-elastic component provided by the model. The results for the three cases were presented in Figures 10-12, respectively (details were mentioned within the figure captions).
thus that the dynamic behavior cannot be differently influenced by the fractional derivative than the Voigt-Kelvin-based model. Hereby, in the frame of this study and taking into account the nonlinear damping expression, becomes interesting to analyze a VdP-type model, embedding fractional order derivative, related to its capability that emulate a classical VdP model. The following expression of VdP nonlinear model including α-fractional derivative was proposed for computational investigations [30]   where f(t) denotes the response of the system, x(t) is the instantaneously displacement, and c0, k0 denote damping and stiffness coefficients, respectively. In respect of the hypothesis within previous case analysis, the parameters c0, k0 were identically valued as follows: c0 = 1, k0 = 0. Taking into account the nonlinear term in parenthesis within the left side hand of Equation (22), the analysis has to be conducted related to displacement magnitude. In this study, random values for magnitude were taken. The results presented correspond to the 3.3, 1.0, and 0.3 valued displacement magnitudes. The monitored parameters, for each case, were similarly adopted than the previous linear case: Timed evolution, spectral composition, hysteretic loops with Poincare maps overlapped, and hysteresis axis according to the pseudo-elastic component provided by the model. The results for the three cases were presented in Figures 10-12, respectively (details were mentioned within the figure captions).  A comparative analysis between this group of results highlight the fact that the VdP model embedding the fractional derivative is also able to emulate a pseudo-conservative behavior, depending on the α-order. This aspect is dependent to the excitation magnitude and appears for any of its values except the unitary ones-see Figure 11(c,d).
As was performed for the previous linear case, the authors developed a computational application in order to identify the proper formulation of the equivalent VdP visco-elastic model and to evaluate the parameters values of equivalent models. Identification procedure was imposed by the nonlinear formulation of the dissipative term in VdP model and by the fact that this expression has propagated to the conservative term for the equivalent model. In this view, the authors A comparative analysis between this group of results highlight the fact that the VdP model embedding the fractional derivative is also able to emulate a pseudo-conservative behavior, depending on the α-order. This aspect is dependent to the excitation magnitude and appears for any of its values except the unitary ones-see Figure 11c,d.
As was performed for the previous linear case, the authors developed a computational application in order to identify the proper formulation of the equivalent VdP visco-elastic model and to evaluate the parameters values of equivalent models. Identification procedure was imposed by the nonlinear formulation of the dissipative term in VdP model and by the fact that this expression has propagated to the conservative term for the equivalent model. In this view, the authors supposed two forms of the equivalent model. The first case was in respect of the classical VdP expression as follows in the same time that the second case supposes a specific nonlinear terms for equivalent rigidity as follows The analyses were performed for different values of the α-fractional order and displacement magnitudes, respectively. For presentation, the 1.8, 1.0, and 0.8 valued displacement magnitudes were adopted, maintaining α = 0.7. The results were depicted as pair diagrams of the first and second equivalent model, for each displacement magnitude, in Figures 13-15, respectively.
Comparative analysis of this last group of results (depicted in Figures 10-12) show that identification of the right equivalent model formulation is easily done, taking into account the differences between the loops evolutions in diagrams. Even if the first equivalent model (emulated VdP-based model with linear rigidity) is able to provide the same energy dissipation for a proper damping coefficient and supplies certain pseudo-stiffness, the hysteretic evolution is clearly different. This fact can produce certain difficulties on the tuning process between theoretical behavior and the experimental investigations.
In the same time, the second approach (emulated VdP-based model with propagated nonlinear rigidity term) provides the best approximation with the α-fractional derivative model. Obviously, similar to the case of the linear model, the estimated relative error tends to zero, only and directly depending by the computation precision. supposed two forms of the equivalent model. The first case was in respect of the classical VdP expression as follows in the same time that the second case supposes a specific nonlinear terms for equivalent rigidity as follows The analyses were performed for different values of the α-fractional order and displacement magnitudes, respectively. For presentation, the 1.8, 1.0, and 0.8 valued displacement magnitudes were adopted, maintaining α = 0.7. The results were depicted as pair diagrams of the first and second equivalent model, for each displacement magnitude, in Figures 13-15, respectively.
Comparative analysis of this last group of results (depicted in Figures 10-12) show that identification of the right equivalent model formulation is easily done, taking into account the differences between the loops evolutions in diagrams. Even if the first equivalent model (emulated VdP-based model with linear rigidity) is able to provide the same energy dissipation for a proper damping coefficient and supplies certain pseudo-stiffness, the hysteretic evolution is clearly different. This fact can produce certain difficulties on the tuning process between theoretical behavior and the experimental investigations.
In the same time, the second approach (emulated VdP-based model with propagated nonlinear rigidity term) provides the best approximation with the α-fractional derivative model. Obviously, similar to the case of the linear model, the estimated relative error tends to zero, only and directly depending by the computation precision. supposed two forms of the equivalent model. The first case was in respect of the classical VdP expression as follows in the same time that the second case supposes a specific nonlinear terms for equivalent rigidity as follows The analyses were performed for different values of the α-fractional order and displacement magnitudes, respectively. For presentation, the 1.8, 1.0, and 0.8 valued displacement magnitudes were adopted, maintaining α = 0.7. The results were depicted as pair diagrams of the first and second equivalent model, for each displacement magnitude, in Figures 13-15, respectively.
Comparative analysis of this last group of results (depicted in Figures 10-12) show that identification of the right equivalent model formulation is easily done, taking into account the differences between the loops evolutions in diagrams. Even if the first equivalent model (emulated VdP-based model with linear rigidity) is able to provide the same energy dissipation for a proper damping coefficient and supplies certain pseudo-stiffness, the hysteretic evolution is clearly different. This fact can produce certain difficulties on the tuning process between theoretical behavior and the experimental investigations.
In the same time, the second approach (emulated VdP-based model with propagated nonlinear rigidity term) provides the best approximation with the α-fractional derivative model. Obviously, similar to the case of the linear model, the estimated relative error tends to zero, only and directly depending by the computation precision.

Discussion
Each case of analysis was briefly discussed on the same time with its presentation. Besides this, the main findings within this study has to be detailed, in a general view, taking into account the intrinsic availability of α-fractional derivatives that virtually modify the behavioral characteristics of linear/nonlinear vibration isolators models. The results in Figure 8, in terms of equivalent damping and stiffness, respectively, are relevant into this way. From the diagram of damping evolution in respect of the α-order results, the fractional derivative model is able to provide a full scale of dissipation, from weak to strong, related to the full range of α. In the same time, from the diagram of pseudo-elastic component results that also provided a relative very large scale for rigidity, from stiff, by very stiff (more of doubled initial value), to zero stiffness. These two observations are very useful, e.g., for investigations of elastomers-based vibration isolators, which provide strong nonlinear evolutions under the external dynamic loads and involve many parameters within computational (experimentally augmented) models.
Taking into account the comparative analysis of diagrams in Figure 8, it has to be mentioned that the evolutions of the two components are strongly correlated. In this sense, for small α orders, the FDE-based model enables only weak-dissipation with stiff to very stiff behavior. For medium to 0.9 values of α, the FDE-based model still enables a stiff behavior with damping transition from medium to strong. Finally, for 0.9...1.0 valued α, the model provides only a soft to zero stiffness, but very damping, behavior. These last observations have restricted the area of practical using for the FDE-based dynamic models. Nevertheless, these theoretical models are able to provide enough of a range of dynamic characteristics, in respect of a single parameter, and was thus suitable for inclusion into the computational approaches of vibration isolator dynamics.
In addition, the availability of the α-fractional derivatives has to be underlined, embedded to the differential constitutive equations, which do not change the hysteretic evolution of a dynamic system with symmetrical evolution in respect to its static state equilibrium point. This is a very useful observation, because the investigations done by the authors highlights that the insertion of fractional-order derivatives into the classical differential formulations of dynamic systems actually provides the preservation of the symmetry of the hysteretic behavior.
For the simple, rubber-like, composites without any augmentations related to differentiation of their reversible evolution under dynamic loads, the previously mentioned fact focuses on the interest to these computational approaches, because of theirs accurate capability to simulate the

Discussion
Each case of analysis was briefly discussed on the same time with its presentation. Besides this, the main findings within this study has to be detailed, in a general view, taking into account the intrinsic availability of α-fractional derivatives that virtually modify the behavioral characteristics of linear/nonlinear vibration isolators models. The results in Figure 8, in terms of equivalent damping and stiffness, respectively, are relevant into this way. From the diagram of damping evolution in respect of the α-order results, the fractional derivative model is able to provide a full scale of dissipation, from weak to strong, related to the full range of α. In the same time, from the diagram of pseudo-elastic component results that also provided a relative very large scale for rigidity, from stiff, by very stiff (more of doubled initial value), to zero stiffness. These two observations are very useful, e.g., for investigations of elastomers-based vibration isolators, which provide strong nonlinear evolutions under the external dynamic loads and involve many parameters within computational (experimentally augmented) models.
Taking into account the comparative analysis of diagrams in Figure 8, it has to be mentioned that the evolutions of the two components are strongly correlated. In this sense, for small α orders, the FDE-based model enables only weak-dissipation with stiff to very stiff behavior. For medium to 0.9 values of α, the FDE-based model still enables a stiff behavior with damping transition from medium to strong. Finally, for 0.9...1.0 valued α, the model provides only a soft to zero stiffness, but very damping, behavior. These last observations have restricted the area of practical using for the FDE-based dynamic models. Nevertheless, these theoretical models are able to provide enough of a range of dynamic characteristics, in respect of a single parameter, and was thus suitable for inclusion into the computational approaches of vibration isolator dynamics.
In addition, the availability of the α-fractional derivatives has to be underlined, embedded to the differential constitutive equations, which do not change the hysteretic evolution of a dynamic system with symmetrical evolution in respect to its static state equilibrium point. This is a very useful observation, because the investigations done by the authors highlights that the insertion of fractional-order derivatives into the classical differential formulations of dynamic systems actually provides the preservation of the symmetry of the hysteretic behavior.
For the simple, rubber-like, composites without any augmentations related to differentiation of their reversible evolution under dynamic loads, the previously mentioned fact focuses on the interest to these computational approaches, because of theirs accurate capability to simulate the realistic behavior, experimentally acquired, in the same time with keeping the nonlinear characteristics identified to the microstructural and phenomenological analyses.
For the sake of comparison in the view of previous paragraph remarks, the diagrams in Figure 16 show the evolution of a rational-power simulated system, with quite identical material parameters to those used within this study. It was adopted in the force to displacement graphs, considering both the conservative and dissipative force components, in order to highlight the unsymmetrical hysteresis of such computational models. The conservative component of the model was supposed to have linear expression to displacement, with a constant value of rigidity. realistic behavior, experimentally acquired, in the same time with keeping the nonlinear characteristics identified to the microstructural and phenomenological analyses.
For the sake of comparison in the view of previous paragraph remarks, the diagrams in Figure  16 show the evolution of a rational-power simulated system, with quite identical material parameters to those used within this study. It was adopted in the force to displacement graphs, considering both the conservative and dissipative force components, in order to highlight the unsymmetrical hysteresis of such computational models. The conservative component of the model was supposed to have linear expression to displacement, with a constant value of rigidity. 1.7] (c). It was assumed a single time-period of 1 Hz excitation signal, and ratio between rigidity to damping, for linear components, valued to 20 s -1 . Within these diagrams, the symbol α must be considered as the rational power of additional dissipative component (multiplied by 20 from linear damping coefficient) into the computational model, with specific values mentioned on graph legends.

Conclusions
The computational investigations within this study highlight a very interesting aspect of the α-fractional derivatives embedded in classical linear or nonlinear oscillators. This aspect consists by the capability of emulating a pseudo-conservative characteristic for a computational implementation of the only dissipative component. In the view of practical using these models for the evaluation of vibration isolators dynamic behavior, this capability is very useful because it provides a single tuning parameter (in term of the α-order) instead of two (for linear case) or more parameters (for nonlinear case), for classical rheological dynamic models. Moreover, these FDE-based approaches are able to maintain the hysteresis symmetry, specific to the real systems, within its dynamic behavior under external loads.
Two cases of classical rheological models were analyzed in this work, fractional-derivative augmented. The future investigations will be focused on enlarging the area of fractional calculus applicability within other basic rheological models, and the potential useful results will be discussed into forthcoming works.
Author Contributions: The all authors contributed equally and all authors read the manuscript and approved the final submission.
Funding: This research received no external funding.

Conclusions
The computational investigations within this study highlight a very interesting aspect of the α-fractional derivatives embedded in classical linear or nonlinear oscillators. This aspect consists by the capability of emulating a pseudo-conservative characteristic for a computational implementation of the only dissipative component. In the view of practical using these models for the evaluation of vibration isolators dynamic behavior, this capability is very useful because it provides a single tuning parameter (in term of the α-order) instead of two (for linear case) or more parameters (for nonlinear case), for classical rheological dynamic models. Moreover, these FDE-based approaches are able to maintain the hysteresis symmetry, specific to the real systems, within its dynamic behavior under external loads.
Two cases of classical rheological models were analyzed in this work, fractional-derivative augmented. The future investigations will be focused on enlarging the area of fractional calculus applicability within other basic rheological models, and the potential useful results will be discussed into forthcoming works.
Author Contributions: The all authors contributed equally and all authors read the manuscript and approved the final submission.
Funding: This research received no external funding.