#### Abstract

Viscoelastic materials are widely used in structural dynamics for the control of the vibrations and energy dissipation. They are characterized by damping forces that depend on the history of the velocity response via hereditary functions involved in convolution integrals, leading to a frequency-dependent damping matrix. In this paper, one-dimensional beam structures with viscoelastic materials based on fractional derivatives are considered. In this work, the construction of a new equivalent viscous system with fictitious parameters but capable of reproducing the response of the viscoelastic original one with acceptable accuracy is proposed. This allows us to take advantage of the well-known available numerical tools for viscous systems and use them to find response of viscoelastic structures. The process requires the numerical computation of complex frequencies. The new fictitious viscous parameters are found to be matching the information provided by the frequency response functions. New mass, damping, and stiffness matrices are found, which in addition have the property of proportionality, so they become diagonal in the modal space. The theoretical results are contrasted with two different numerical examples.

#### 1. Introduction

The relevance of the correct modeling of damping mechanisms in structural dynamics and in particular mechanical vibrations is well known. Directly related to this point, the linear viscoelastic models have been playing a leading role in the last decades. On one hand, the linearity allows us to use very efficient tools for the numerical solution. On the other hand, the great variety of mathematical models covers almost the totality of the materials and structures used in engineering and scientific applications within vibrations and control. In the last decades, a great number of authors focused their research on enhancing the mathematical models to solve new technological challenges in aerospace, civil, and mechanical engineering, for instance, new vibration isolation and control systems or the use of frequency-dependent materials. Thus, the works of Flugge [1], Nashif [2], Jones [3], and more recently Adhikari [4, 4] have been of special importance, including the main contributions to this science. In general, a viscoelastic model is characterized by dissipative forces depending on the time history of the response through a convolution integral involving kernel hereditary functions. Therefore, the equations of motion for a multiple degrees of freedom system is presented as an integrodifferential equation, saywhere are the mass and elastic stiffness matrices, respectively. For continuous systems, the Finite Element Method will provide the corresponding tools to determine such matrices. On the other hand, contains the viscoelastic kernel functions in time domain and represents the array (in column) with the degrees of freedom. Each viscoelastic model is characterized by a kernel function; a survey of such models can be found in [5, 6]. Due to the linearity, the response of these models has a straightforward representation in frequency domainwhere and are the response and excitation Fourier transforms, respectively. contains the Fourier transform of the kernel functions. However, the time-domain response becomes very expensive computationally, especially for large systems. Several works have been focused on its resolution; among them, those ones based on the state-space approach result in special interest. Thus, the works of Golla and Hughes [7] and McTavish and Hughes [8] introduced the Golla-Hughes-McTavish (GHM) method based on using new internal variables. Muravyov [9, 10] and Muravyov and Hutton [11, 12] formulated a general Laplace domain-based method for isotropic and homogeneous hereditary materials with exponential kernels, analyzing both the forced and the free vibration response. Menon and Tang [13] proposed a solution for viscoelastic problems with multiexponential kernels based on transforming the characteristic equation into a polynomial. Thus, the order of the system was raised as many times as the number of exponential kernels.

Among the mathematical models proposed in viscoelasticity, those that are based on the fractional calculus are of special significance for this paper. In 1983, Bagley and Torvik [14, 15] presented the theoretical foundations for the viscoelastic models based on fractional derivatives. From the time-domain point of view, the theory of hereditary solid materials was generalized by Koeller [16] to the fractional calculus, showing this application to different spring-damper systems. Gaul et al. [17] have used the fractional calculus in the constitutive equations, describing a particular case of analytical integration. They compared the results with different nonfractional models with the objective of obtaining equivalent damping ratios. Gaul [18] studied the influence of the viscoelastic models based on fractional and integer derivatives in wave propagation and structural vibration, analyzing the stress relaxation and the creep of strain. Pritz [19] analyzed in depth the four-parameter model for real solid materials, describing the influence of the loss factor peak in the viscoelastic behavior. In subsequent works of this author [20–22], the frequency dependence of the different constitutive parameter is examined with detail, especially its application to rubbers, rubber-like materials, and elastomers. Rossikhin and Shitikova [23] published a survey of the most important advances in the application of fractional derivatives applied to dynamical problems.

Our main motivation is to transform the system governed by (1) into a simpler one, much more accessible from a numerical point of view. Traditionally, the most popular dynamical system involves viscously damped models, whose dissipative forces are proportional to the velocity of the degrees of freedom. Thus, the challenge is to find an equivalent viscous model, for the particular case of one-dimensional viscoelastic structures, defining not only an equivalent damping matrix but also a new mass matrix and a new stiffness matrix , so that the response given by the equationscan be considered as a very reasonable approximation to the real one, corresponding to the viscoelastic system. The search of equivalent dynamic models capable of representing the response of the original one with enough accuracy is common in structural dynamics. Bandstra [24] obtained for single-degree-of-freedom systems (SDOF) the equivalent viscous damping for different nonlinear systems with a method based on the equivalence of the energy dissipation. For one-dimensional structures, Banks [25] used the modal analysis for extracting the set of the modal damping ratios. These ones were deduced for the special case of viscous-air damping combined with a constitutive relationship based on a Kelvin-Voigt model. Concerning more complex structures with added viscoelastic dampers, various methods have been proposed for their modeling through equivalent viscous models. Lima [26] suggested an efficient methodology for assembling structural systems with viscoelastic dampers based on a frequency response function coupling technique. Ungar and Kerwin [27] proposed in 1962 the Modal Strain Energy method to find the modal damping ratio in a general viscoelastic system. This method was implemented in a Finite Element Program by Johnson and Kienholz [28]. Bilbao et al. [29] proposed a proportional damping matrix, using optimization based methods. Genta and Amati [30] studied the hysteretic models in rotor dynamics and their equivalent damping ratios, analyzing the advantages and the disadvantages of using them for obtaining the time-domain response. In the specific case of a viscoelastic system based on fractional calculus, Makris [31] proposed an equivalent viscous model giving an equivalent damping ratio and a new natural frequency for viscoelastic dampers in single-degree-of-freedom systems. The theoretical problem of finding a second-order model equivalent to a general viscoelastic system has been studied by Segalman [32], assuming small viscoelasticity. In this work, an equivalent second-order model was deduced, finding a damping matrix and a new stiffness matrix. However, the resulting model is not viscous due to the involved matrices being in general defined in the complex domain. Adhikari [33] has proposed also a second-order system in an interesting work, where the key idea is to evaluate the viscoelastic function in each complex frequency, assuming also low viscoelasticity. For extracting the set of complex eigenvalues, Adhikari and Pascual [34] proposed an efficient algorithm based on the expansion of the viscoelastic function.

In this paper, a new viscous approach valid for isotropic one-dimensional viscoelastic structures based on the fractional derivative is proposed. The method allows us to find three system matrices forming a new equivalent viscous system (3), say mass, damping, and stiffness, capable of representing the response of the original viscoelastic model given by (1). Two main ideas have been highlighted in the construction of this model: first, viscous damping ratio and natural frequency arise after imposing equality between complex eigensolutions; second, a new equivalent mass can be found by forcing magnitude of the frequency response function to be equal around the damped frequency. The proportional nature of such systems is precisely the key condition to derive the method. Obviously, the available tools to solve viscous system have in general many more advantages, both in computational efficiency and in mathematical simplicity. Furthermore, the obtained matrix is proportional, that is, becomes diagonal in the modal space of the mass and the stiffness matrices.

As viewed in the bibliography review, works concerning the proposal of equivalent models for viscoelastic systems are in general only addressed to find equivalent damping matrices. Modifications in the dynamics matrices, say mass and stiffness, are usually not introduced, especially in the mass matrix. In addition, lightly viscoelasticity becomes a commonly assumed hypothesis in the majority of the cases. The motivation of this paper is double: not only to build a proportional viscous model whose response can be a good approximation to the exact one, given by the viscoelastic model but also extending the resulting model to systems with a relatively high viscoelasticity.

#### 2. Viscoelastic Constitutive Relationships

In this point, the mathematical model of viscoelastic damping forces will be presented. As stated in Introduction, one-dimensional nonviscously damped structures are considered. The general forms of constitutive relationships for these materials based on fractional derivatives were introduced by Bagley and Torvik [15] and they consist in linear fractional differential equations relating stresses and strains. The same authors suggested in [35] a simpler relation, with relatively few terms allowing describing a wide range of materials. They also demonstrated that the order of the fractional derivatives involved in the equation must be the same for satisfying the thermodynamic constraints. The resulting model is named* four-parameter model* and it has been used in this paper for modeling the constitutive relationships in one-dimensional structures, which have the following form:where represent plane stresses and strains in any point of the cross section of the beam. The operator denotes the fractional derivatives of order . In the application of the fractional derivatives to viscoelastic models, it is assumed that [14]. The parameter with dimensions of time is named* relaxation time*. and are, respectively, the static Young’s modulus and shear modulus, that is, the dynamics modulus at the null frequency (). Parameters and are the high frequency limit value of the dynamic modulus (), also named* storage modulus* because in general it is verified that and with . The storage modulus is equal for both longitudinal and shear behavior, something that is compatible with the fact that Poisson’s ratio presents in general small variations with frequency for real materials [36–38]. Pritz [19] studied the four-parameter model describing the influence of the parameters in the frequency domain for a great variety of real solids materials.

Using the properties of the fractional derivatives, (4) can be transformed into the frequency domain. The differential relations become now algebraic. Thus,where and are the complex Young modulus and the complex shear modulus, respectively. After some straight operations, their expressions result inTherefore, only three parameters, say , , and , will be required to describe the behavior of the dynamic modulus in frequency domain through the following complex function:As it will be seen later, it sometimes becomes interesting to know if our parameters induce low or high viscoelasticity. Thus both the real and the imaginary parts of the functions and provide valuable information. The frequency dependence in both complex moduli is known as* dispersion* and is controlled by the function . Solving for the real and imaginary parts, this function can be expressed bywhere the functions and denote the dynamic and loss modulus, respectively. The viscoelasticity of the model is directly related with the variation of in the frequency domain. Strong variation is associated with high viscoelasticity, while low viscoelasticity is a property of materials whose complex modulus does not present great dependence of the frequency. The function is the named* loss factor* of the material. This function can be used for measuring the level of viscoelasticity of the system [39]. Pritz [19, 40] studied this function for real solid materials providing two interesting properties. First, the loss factor of a three-parameter model like (7) presents one peak, in agreement with the behavior of the majority of materials. Second, the value of the loss factor at the peak , named* loss factor peak*, is a good indicator of the quantification of the viscoelasticity in the dynamic model. The expression of the loss factor peak as a function of the parameters of the damping model isAccording to Pritz, a low loss factor peak, , is characteristic of hard plastics and other structural materials (wood, concrete, metals, etc.) resulting in the weak frequency dependence of the dynamic modulus. On the other hand, the large loss, , is characteristic of rubbers and rubber-like materials used for vibration isolation and damping. The loss factor peak will be used later in the numerical examples due to its relevance in the results. Other indicators of the quantification of the viscoelasticity are available in the work of Adhikari and Woodhouse [41].

#### 3. Equations of Motion in the Frequency Domain

The governing equations of one-dimensional structures with constitutive relationships based on (4) have the form of a system of fractional differential equations. Although several methods exist in the literature [23], both numerical and analytical, to solve these equations in time domain is far from this scope. Otherwise, the proposed approach in this paper is based on the construction of an equivalent frequency response function. Hence, it becomes necessary to express the motion equations in the frequency domain. For that reason, the principle of virtual work in its frequency domain version [42, 43] will be used. Thus, we can writewhere and represent stress and strain arrays in frequency domain, is the density, and is the circular frequency (in rad/s). The vector is the Fourier transform of the displacements field in the volume , where represents the interval where the beam is defined and is the cross section. contains the spectral displacements and rotations of the axis. The vectors and represent the Fourier transforms of the distributed loading and the applied load vector, respectively. The longitudinal domain is discretized in finite elements. The solution over the longitudinal -axis can be expressed as function of the variables in the nodes bywhere denotes the matrix containing the shape functions. The kinematics equations relate displacements with the strains field and with the vector through certain matrices defined for each point on the cross section [44].whereIntroducing the relations (13) in the principle of virtual work (11) and using the constitutive relationships in the frequency domain result in the following system of equations:whereare the mass and the stiffness matrices of the structure and is the size of the vector . From (15), the frequency domain response can be found asresulting in the receptance matrix . The next section describes the general methodology of the proposed method in this paper. Since the complex eigenvalues of the viscoelastic system (15) play an important role in this development, their numerical computation is analyzed as a separated point.

#### 4. Proposal of the Viscous Model

##### 4.1. Complex Eigenvalues

Let us consider the eigenvalue problem corresponding to the homogeneous system in (15).In general, the number of eigenvalues of this problem can be greater than , say , depending on the nature of the viscoelastic model [45]. The additional roots are named* nonviscous modes* and are in general overcritically damped. In this paper, we only consider the * elastic modes* corresponding to the complex roots. Since the set of complex eigenvalues is formed by complex-conjugate pairs in the Laplace domain, the damped eigenfrequencies can be ordered in a set with the form , where denotes the complex-conjugate. Obtaining the set of complex eigenvalues of problem (15) is a key condition in the developing of the method.

As a previous step, it is necessary to solve the undamped elastic eigenvalue problem. Therefore, let us consider by the th mass-normalized th eigenvector and is its undamped natural frequency. Arranging eigensolutions in the matrices and , both within , the following relations are well known from the classical modal analysis:where is the th-order identity matrix. The change of variable transforms in a diagonal system the eigenvalue problem of (18). Hence, the complex eigenvalues are the solution of the equationIn particular, the th complex eigenfrequency must verify the equationwhere is the th undamped natural frequency. Adhikari and Pascual [34] studied the numerical solution of proportional or lightly nonproportional viscoelastic systems, proposing a method for the following equation, in the Laplace domain ():The key idea is the expansion of around , where . It can be easily shown that our equation (21) can be expressed as the form (22), taking a viscoelastic function defined bywhere is the function defined by (7) but changing . However, in this particular case, a singularity appears when the previous limit is evaluated. Actually, applying L'Hôpital's rule results inassuming that and . Therefore, it is not suitable to apply the method just as it is described in [34]. However, the expansion of is always possible provided that another valid initial point can be found. In this work, the following alternative initial point is proposed:This approximation comes from solving the characteristic equation after evaluating the viscoelastic function at the th undamped natural frequency; that is,The negative root has not been considered, since only the complex solution with a positive real part is of interest in this step. In general, it is assumed that the complex function presents a smooth variation in frequency domain [20] and light variations around are expected as well. As a consequence of that, expanding in its Taylor series up to the second order should be a good approximation for the complex frequencies. Thus, if is a root, there exists a number , such that . Thus, a second-order polynomial results:with the following complex roots:whereAmong both roots of (28), that one with positive real part is chosen, assuming that such value will be closer to . This numerical approach will be used in the rest of this paper to find eigenvalues of the viscoelastic systems.

##### 4.2. Single-Degree-of-Freedom Systems

For the sake of clarity, it is interesting to present first the fundamentals of the proposed method for single degree of freedom (SDOF). Figure 1 shows schematically a SDOF system with a mass attached to a fixed point through a viscoelastic constraint modeled with fractional derivatives.

Let and denote the force reaction and the time-dependent external force, respectively. The relationship between and the degree of freedom can be expressed as (see (4))where , , and are the parameters of the viscoelastic model based on fractional derivatives. Applying the Fourier transform to this expression, it follows that , where denotes the complex stiffness [2]. Assuming an excitation force , the motion equation in frequency domain can be expressed asand we introduce as the frequency response function (FRF) of the system:where is the undamped natural frequency of the system. Now, using as an initial point(28) can be used to find a closed-form approximation of the eigenvalue . Of course, represents an approximation of the other root (associated with conjugate-complex pair).

The aim is to build a linear and strictly viscous model using for that the information provided by the viscoelastic original one. Therefore, only the three parameters, say mass , viscous damping , and stiffness , are required for its definition. As is well known, the FRF of that model will have the following form:where the* equivalent undamped natural frequency* and its* equivalent damping ratio* have been introduced. Our first condition is to impose the fact that the poles of are the same as those of the original FRF, . Consequently, the main feature of the approached function is to take the form of a second-order polynomial with roots and . Thus,Now, the complex frequencies will be expressed as the product of a nondimensional complex number and the undamped natural frequency associated. Hence, the following can be written:where are two nondimensional parameters. It is worth highlighting that is analytically available from (28) and therefore both the real part and the imaginary part are closed-form expressions explicitly defined as functions of the parameters of the original model. Mathematically, they can be expressed as certain functions and such thatIntroducing (36) into (35) and expanding the product result inIdentifying the last two expressions, the damping factor and the natural frequency of the equivalent viscous model can be obtained as functions of and .Secondly, the proposed viscous approach is completed calculating an equivalent mass using information provided by the frequency response functions. The key idea is to compare the frequency response functions around the peaks. As is well known, for a viscous model with damping factor and natural frequency , the peak of the absolute value of the FRF is located in the abscissa . Due to the relevance of the range around this point in the response, it is reasonable to calculate the mass under the condition that both FRFs have the same value at . Mathematically, this can be expressed by the equality between both models.Therefore, the mass of the equivalent model can be expressed byUsing (39), results:Substituting this expression in and and after some straight operations, a more compact expression for can be deduced:To sum up, it has been proven that the three equivalent parameters, , are directly related to the parameters of the original model , and through the relationsand using (37). Above, the notation has been introduced. In this point, essentials of the proposed method have been presented using SDOF systems. The next point is focused on the generalization of these concepts to multiple-degrees-of-freedom systems (MDOF).

##### 4.3. Multiple-Degrees-of-Freedom Systems

As is well known, the entries of receptance matrix of a multiple-degrees-of-freedom (MDOF) system are the FRFs. As shown in (15), the equations of motion in frequency domain can be deduced from the inversion of the coefficients matrix for each .In general, to calculate as many inverse matrices as the number of frequencies is very expensive computationally. Very efficient tools are available in the modal analysis to find the receptance matrix [46]. They are based on the modal space of the generalized eigenvalue problem of and . Hence, a matrix exists, such thatwhere and are two diagonal matrices whose entries are related to the natural frequencies through the equation . The mass-normalization of each eigenvector shows as a result the matrix previously presented in (19). Thus, the receptance matrix is obtained from the following matrix product:Each entry of the main diagonal in is the FRF of the th modal coordinate, . The resulting FRF has the same form as that of a SDOF system defined in (32) with mass and stiffness . Thus, introducing again the function , which only depends on the parameters of the damping model and on the th undamped frequency, the th modal FRF leads toThe set of complex frequencies are computed using the above approach (28). Thus, as a simple extension of the method described for the SDOF systems, the function can be approximated bywhere the equivalent parameters , , and are given bywith and . These equations can be put in matrix form bywhere and and the two diagonal matrices, and , have been introduced, which are working as transformation matrices in the modal space:Therefore, the approximation of is given byFurthermore, the receptance matrix of the MDOF system can be estimated with

The above expression is just the classical form of a receptance matrix corresponding to a linear-viscous model. Now, three matrices , , and can be defined so that they verify the following conditions:Finally, the proposed viscous approach for MDOF viscoelastic systems is defined as the linear system and has of (55) as receptance matrix. Immediately, the relations (56) show that the mass, damping, and stiffness matrices are diagonal in the same base. Moreover, the three matrices are symmetric, and therefore Maxwell’s Reciprocity Theorem is also verified in the proposed model. After some straight algebra, one finds closed-form expressions for and . Thus, starting from their definition and using (46) and (52), we havewhere denotes the inverse transpose of . Above, the following matrices have been introduced:which can be interpreted as transformation matrices affecting the original matrices and to provide as an output their equivalent matrices and , respectively.

The new viscous model obtained from (58) is considered as the main contribution of this paper. It should be noted that the new matrices are somehow fictitious; they have been found as result of a mathematical manipulation. However, it is expected that this new fictitious viscous model will estimate accurately the response of the original nonviscous model provided that its receptance matrix, , is a good approximation of the original model one, . It must be highlighted that the derivations have been carried out in the frequency domain and the method has been built mainly trying to match the transfer functions of viscoelastic and proposed models. As consequence, good behavior of steady-state response of harmonically forced systems is expected when initial transient time range is negligible. Thus, problems in which transient response is relevant (for instance, seismic analysis) fall out of the scope of the proposed methodology.

To illustrate the efficiency of the method, a set of different numerical examples will be examined in the next point. It will be shown that the limitations of the proposed method are directly related to the viscoelasticity of the system, measured with the loss factor given in (10).

#### 5. Numerical Examples

##### 5.1. Example 1: Single-Degree-of-Freedom Systems

In this point, the application of the theoretical results seen above for SDOF will be investigated. Figure 1 shows schematically a SDOF model with a mass kg attached to a fixed point through a linear spring with stiffness kN/m. The viscoelastic model directly depends on the parameters , , and and, indirectly, on the loss factor peak (see (10)). This latter can be considered as a good index to measure the model viscoelasticity. The results are ordered in four different numerical cases (see Table 1), depending on this parameter: Case 1 corresponds to the lowest value of the loss factor peak (low viscoelasticity) and Case 4 to highest one (high viscoelasticity). Additionally, also in Table 1, the results have been compared with those of the equivalent parameters of the proposed method, calculated from (39) and (44).

The results presented in Table 1 show that the higher the viscoelasticity of the system is, the larger the differences between the equivalent parameters and the real ones, and , are. Therefore, it is expected that the loss of accuracy between the responses of the proposed method and the exact one will be highly correlated to the values of the loss factor peak, .

Figures 2(a)–2(d) show the frequency response functions (FRF) in both magnitude and phase corresponding to the four different cases listed in Table 1. The viscoelastic model affects the frequency response moving rightwards the location of the resonance peaks [19, 40]. From Case 1 to Case 4 and as the system becomes more damped, the loss factor peak somehow shifts the FRF curves to the right in the range of frequencies. Our proposed viscous model is able to reproduce this behavior, moving also their FRFs to the right as the damping conditions change. This feature arises from the fact that the equivalent parameters, and , have been calculated ensuring equality between complex frequencies (see (51)). The effect of the equivalent mass is also presented in Figures 2(a)–2(d). Indeed, since is calculated imposing the same magnitude of FRFs (for both viscoelastic and viscous models), it is expected that around the resonance peaks the proposed method estimates accurately the value of the magnitude. Further than the resonance peak, we observe that results for low frequencies range differ from the exact ones as the system becomes more damped. This lack of accuracy can be explained after evaluating and comparing the FRFs at the static problem, . ThusAccording to Table 1, the value increases with . Additionally, since the phase of the response has not been taken into account to form the equivalent viscous model, this variable of the response (phase) presents a lack of precision as the viscoelasticity increases.

**(a)**CASE 1,

**(b)**CASE 2,

**(c)**CASE 3,

**(d)**CASE 4,In order to study the evolution and the range of the accuracy of the proposed viscous model, the relative error between exact and proposed FRFs (in percentage) defined ashas been plotted for the four cases of study in Figure 3. The curves show a clear correlation between the loss factor peak and the error. For low frequencies, every case shows the maximum errors, a fact that has been explained above. As the frequency increases, the benefits of imposing the same values for the FRF magnitude explain a decrease of the error up to a stable value (horizontal asymptotes) that varies between 0.2% and 10%, depending on the damping case. The main source of error comes from the differences shown in phase. These discrepancies, although bounded, seem to remain constant along the frequency range, as shown in Figure 2 (Case 4). On the other hand, the magnitude of the FRF (in the denominator of (61)) tends to decrease with frequency. Thus, it is expected that the error does not decrease with frequency, something that can be observed in Figure 3.

##### 5.2. Example 2: Viscoelastic 1D Structures

This section is focused on the validation of the proposed method for fractional derivative-based viscoelastic damping 1D structures. The beam model described in Sections 2 and 3 is used. The finite element model is formed by finite elements. Two degrees of freedom are associated with each node: the vertical displacement and the rotation . Boundary conditions of the type clamped-pinned are considered, as shown in Figure 4. The mechanical properties considered for the numerical example are summarized in Table 2.

As for the first example, the viscoelastic parameters will be presented and ordered in Table 3, considering the four numerical cases, according to the level of viscoelasticity, measured by the loss factor peak, . As a first step, the set of complex frequencies , need to be calculated. They will be approximated using the proposed approach given by (28). It is expected that both the real part and the imaginary part show relative errors with respect to the exact eigenvalues (computed using iterative methods [47]). In order to compare and check the accuracy of the eigenvalues approximation, maximum errors within the array containing the entire set of complex frequencies will also be calculated. Let and denote the maximum error in real imaginary part, respectively. Mathematically, they can be found bywhere represents the total number of modes considered. These values are listed in the two last columns of Table 3. Even for highly damped systems, the numerical approach presents very low relative errors, validating the numerical method. As proved, the construction of the equivalent viscous model depends strongly on this approximation so that the better the fit of eigenvalues, the higher the accuracy of the frequency response function.

Frequency response function (FRF) of the exact model, , and its viscous approximation, , can be found using (47) and (56), respectively. The 4 numerical cases of Table 3 have been plotted in Figures 5(a)–5(d). The viscous approach proposed in this paper includes a new mass matrix, arising from the condition of forcing the magnitude of FRFs of both viscoelastic and viscous models to be equal at the resonance frequencies (see (42) and (43)). Therefore, it will be expected that magnitude will present high accuracy along the complete range of frequencies of study. This prediction can be visualized in Figures 5(a)–5(d), where even for highly damped beams, the error in the absolute value of FRF remains very low. It can also be observed that phase plots fit very well with exact curves despite presenting more visible differences than those of the magnitude curves, especially as the viscoelasticity increases (Cases 3 and 4). Finally, a discrepancy in the phase plots at frequencies between the second and the third modes is appreciated in all phase plots of Figures 5(a)–5(d). The explanation of this jump in the phase must be searched in the interpretation of as a curve in the complex plane, controlled by the parameter . When the curve goes through the real axis (negative part), from the third to the second quadrant, a jump in the phase arises from to because of the convention . Considering these comments, the actual error of around this point does not present great differences with respect to those at the proximity, something that can be observed in the error plots described in the next paragraph (see Figure 6).

**(a)**CASE 1,

**(b)**CASE 2,

**(c)**CASE 3,

**(d)**CASE 4,As it has been seen before, it is interesting to compare the four damping cases listed in Table 3 through relative error (%) curves along the frequency range, shown in Figure 6. Certain similarities between the behavior of the error curves for multiple DOF and those of single DOF in Figure 3 can be observed. The first one is the fact of having a higher error in the low frequency range before reaching the first resonance frequency. As for single DOF systems, these discrepancies, which are close to the static problem (), arise from the difference between FRFs at . In fact, at , it yields . In the paper, it is also demonstrated that differences between matrices and become higher as the viscoelasticity increases (from Cases 1 to 4). Figure 3 also depicts that the errors tend to remain stable asymptotically, except at certain points (located at the antiresonances) where abrupt peaks of error are produced. As the system becomes more damped, these peaks tend to be smoothed and the global error increases. Considering any case of study (from Case 1 to Case 4), the orders of magnitude of the FRF error in both single and multiple DOF systems are similar, as simple inspection of Figures 3 and 6.

#### 6. Conclusions

In this paper, one-dimensional beam structures of viscoelastic materials with damping based on fractional derivatives are under consideration. Mechanical properties depend on frequency, leading to a system of integrodifferential equations of motion in time domain. It is well known that viscously damped systems can be solved much more efficiently both in time and in frequency domain. The main motivation in this article is to find an equivalent viscous model that can be able to reproduce numerically the response of the original viscoelastic model with enough accuracy. The developed methodology is based on three stages. Firstly, closed-form expressions, which depend analytically on the viscoelastic damping parameters, are proposed for evaluating the complex eigenfrequencies. Secondly, equivalent viscous parameters (damping ratio and natural frequency) are found under the condition that both models (original-viscoelastic and equivalent-viscous) have the same complex eigenfrequencies. Thirdly, forcing the equality between magnitudes of the frequency response function at the resonance frequencies, a new equivalent modal mass can be found. The approach is described for single-degree-of-freedom systems and extended for 1D beam-like viscoelastic structures with multiple degrees of freedom.

The theoretical developments have been validated using two numerical examples. Example 1 shows how to construct a single-degree-of-freedom viscous system equivalent to a viscoelastic model based on fractional derivatives. The theory predicts that frequency response function magnitude of viscoelastic and viscous model fits accurately, something that can be observed in the numerical results. The level of accuracy strongly depends on the amount of viscoelasticity in the system, measured by the loss peak factor. Example 2 is focused on the validation of the proposed method for 1D viscoelastic beams. Eigenfrequencies have been computed using the proposed numerical approach showing very low relative errors. The complex frequencies prediction affects the quality of the proposed model, in particular the accuracy of the frequency response function. Furthermore, the same behavior as that for single-degree-of-freedom systems is observed: the stronger the damping the less the precision between exact approach and our viscous approach.

#### Data Availability

No data were used to support this study.

#### Conflicts of Interest

The authors declare that there are no conflicts of interest regarding the publication of this paper.