# Proposal of a Viscous Model for Nonviscously Damped Beams Based on Fractional Derivatives

Authors: Mario Lázaro; Jose M. Molines-Cano; Ignacio Ferrer; Vicente Albero
Journal: Shock and Vibration
 (2018)

---

## 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.


---

## Body

## 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, say(1)Mv¨+∫0tGt-τv˙dτ+Kv=Ftwhere M,K∈Rn×n 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, Gt contains the viscoelastic kernel functions in time domain and vt∈Rn 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 domain(2)-ω2M+iωGω+Kv^ω=F^ωwhere v^ and F^ are the response and excitation Fourier transforms, respectively. G(ω) 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 Ce but also a new mass matrix Me≠M and a new stiffness matrix Ke≠K, so that the response given by the equations(3)Mev¨+Cev˙+Kev=Ftcan 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 Ce 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 namedfour-parameter model and it has been used in this paper for modeling the constitutive relationships in one-dimensional structures, which have the following form:(4)σx+Trαdασxdtα=Eϵx+E∞Trαdαϵxdtατxy+Trαdατxydtα=Gγxy+G∞Trαdαγxydtατxz+Trαdατxzdtα=Gγxz+G∞Trαdαγxzdtαwhere σx,τxy,τxz,ϵx,γxy,γxz represent plane stresses and strains in any point of the cross section of the beam. The operator dα/dtα denotes the fractional derivatives of order α. In the application of the fractional derivatives to viscoelastic models, it is assumed that 0<α<1 [14]. The parameter Tr with dimensions of time is namedrelaxation time. E and G are, respectively, the static Young’s modulus and shear modulus, that is, the dynamics modulus at the null frequency (ω→0). Parameters E∞ and G∞ are the high frequency limit value of the dynamic modulus (ω→∞), also namedstorage modulus because in general it is verified that E∞=cE and G∞=cG with c>1. 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,(5)σ^xω=E^ωϵ^xω,τ^xyω=G^ωγ^xyω,τ^xzω=G^ωγ^xzωwhere E^(ω) and G^(ω) are the complex Young modulus and the complex shear modulus, respectively. After some straight operations, their expressions result in(6)E^ω=E1+ciωTrα1+iωTrα≡Eφω,G^ω=G1+ciωTrα1+iωTrα≡GφωTherefore, only three parameters, say α, c, and Tr, will be required to describe the behavior of the dynamic modulus in frequency domain through the following complex function:(7)φω=1+ciωTrα1+iωTrα,φ0=1,(8)φ∞=cAs 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 E^ω=Eφω and G^ω=Gφω provide valuable information. The frequency dependence in both complex moduli is known asdispersion and is controlled by the function φω. Solving for the real and imaginary parts, this function can be expressed by(9)φω=Mdω+Mlω=Mdω1+iηωwhere the functions Md and Ml 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 η(ω)=Ml/Md is the namedloss 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 ηm, namedloss 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 is(10)ηm=c-1sin⁡απ/22c+c+1cos⁡απ/2According to Pritz, a low loss factor peak, 10-3≤ηm≤10-1, 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, 10-1≤ηm≤1, 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 write(11)∫Vs^Tδe^dV-∫Vω2ρd^Tδd^dV=∫0Lq^Tδu^dx+Q^Tδu^Qwhere s^=σ^x,τ^xy,τ^xzT and e^=ϵ^x,γ^xy,γ^xzT represent stress and strain arrays in frequency domain, ρ is the density, and ω is the circular frequency (in rad/s). The vector d^x,y,z,ω∈R3 is the Fourier transform of the displacements field in the volume V=A×0,L, where 0,L represents the interval where the beam is defined and A is the cross section. u^x,ω=u^,v^,w^,θ^y,θ^zT contains the spectral displacements and rotations of the axis. The vectors q^x,ω and Q^ω represent the Fourier transforms of the distributed loading and the applied load vector, respectively. The longitudinal domain [0,L] is discretized in ne finite elements. The solution over the longitudinal x-axis can be expressed as function of the variables in the nodes v^ by(12)u^x,ω≈Nxv^ωwhere Nx∈R5×n denotes the matrix containing the shape functions. The kinematics equations relate displacements with the strains field and with the vector v^ through certain matrices E0y,z,E1y,z∈R3×5 defined for each point (y,z) on the cross section [44].(13)d^x,y,z,ω=E1y,zu^x,ω≈E1y,zNxv^ωe^x,y,z,ω=E0u^+E1du^dx≡By,zN¯xv^ωwhere(14)By,z=E0E1∈R3×10,N¯x=NdNdx∈R10×nIntroducing 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:(15)-ω2M+φωKv^=F^where(16)M=∑e=1ne∫0LeNT∫AE1TρE1dANdx∈Rn×nK=∑e=1ne∫0LeN¯T∫ABTDBdAN¯dx∈Rn×nare the mass and the stiffness matrices of the structure and n is the size of the vector v^. From (15), the frequency domain response can be found as(17)v^ω=-ω2M+φωK-1F^ω≡HωF^ωresulting in the receptance matrix H(ω). 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).(18)-ω2M+φωKv^=0^In general, the number of eigenvalues of this problem can be greater than 2n, say 2n+p, depending on the nature of the viscoelastic model [45]. The p additional roots are namednonviscous modes and are in general overcritically damped. In this paper, we only consider the 2nelastic 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 ω~j,-ω~j∗j=1n, 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ϕj∈Rn the jth mass-normalized jth eigenvector and ωj is its undamped natural frequency. Arranging eigensolutions in the matrices Φ=[ϕ1,…,ϕn] and Λ=diag[ω12,…,ωn2], both within Rn×n, the following relations are well known from the classical modal analysis:(19)ΦTMΦ=In,ΦTKΦ=Λwhere In is the nth-order identity matrix. The change of variable v^=Φy^ transforms in a diagonal system the eigenvalue problem of (18). Hence, the complex eigenvalues are the solution of the equation(20)-ω2In+φωΛ=∏j=1n-ω2+φωωj2≡∏j=1nLjω=0In particular, the jth complex eigenfrequency must verify the equation(21)Ljω=-ω2+φωωj2=0,j=1,…,nwhere ωj is the jth 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 (s=iω):(22)djs=s2+sgs+ωj2=0The key idea is the expansion of dj(s) around s0=-ξωj±ωj1-ξ2, where ξ=lims→0G(s)/2. It can be easily shown that our equation (21) can be expressed as the form (22), taking a viscoelastic function defined by(23)gs=φs-1swhere φs is the function defined by (7) but changing iω→s. However, in this particular case, a singularity appears when the previous limit is evaluated. Actually, applying L'Hôpital's rule results in(24)lims→0⁡gs=lims→0⁡φs-1s=lims→0⁡φ′s=lims→0⁡1sTr1-ααc-11+sTrα=+∞assuming that 0<α<1 and c>1. Therefore, it is not suitable to apply the method just as it is described in [34]. However, the expansion of Ljω is always possible provided that another valid initial point can be found. In this work, the following alternative initial point is proposed:(25)s0=iωjφωjThis approximation comes from solving the characteristic equation after evaluating the viscoelastic function at the jth undamped natural frequency; that is,(26)-ω2+φωjωj2=0⟶ω~j0=+ωjφωj∈CThe 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 ω~j0 are expected as well. As a consequence of that, expanding Lj(ω~j) in its Taylor series up to the second order should be a good approximation for the complex frequencies. Thus, if ω~j is a root, there exists a number ε∈C, such that ω~j=ω~j0+ε. Thus, a second-order polynomial results:(27)Ljω~j0+ε≈Ljω~j0+ε∂Lj∂ωω=ω~j0+12ε2∂2Lj∂ω2ω=ω~j0=0with the following complex roots:(28)ω~j≈ω~j0+2ω~j0-ωj2φ0′±ωj2φ0′-2ω~j02-2ωj2φ0-ω~j0ωj2φ0″-2φ0″ωj2-2where(29)φ0=φω~j0,φ0′=∂φ∂ωω=ω~j0,φ0′′=∂2φ∂ω2ω=ω~j0Among both roots of (28), that one with positive real part is chosen, assuming that such value will be closer to ω~j0. 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). Figure1 shows schematically a SDOF system with a mass m attached to a fixed point through a viscoelastic constraint modeled with fractional derivatives.Figure 1
The original single-degree-of-freedom model with viscoelastic behavior based on the fractional derivatives and its equivalent linear viscous model.Dα(•) denotes the α-order fractional derivative.LetR(t) and F(t) denote the force reaction and the time-dependent external force, respectively. The relationship between R(t) and the degree of freedom u(t) can be expressed as (see (4))(30)R+TrαdαRdtα=ku+cTrαdαudtαwhere α, c, and Tr are the parameters of the viscoelastic model based on fractional derivatives. Applying the Fourier transform to this expression, it follows that R^(ω)=k^(ω)u^(ω), where k^(ω)=φ(ω)k denotes the complex stiffness [2]. Assuming an excitation force F(t), the motion equation in frequency domain can be expressed as(31)-mω2+φωku^=F^ωand we introduce H(ω) as the frequency response function (FRF) of the system:(32)u^=HωF^ω,Hω=1m1-ω2+φωωu2≡1m1Lωwhere ωu=k/m is the undamped natural frequency of the system. Now, using as an initial point(33)ω~d0=ωuφωu=ωu1+ciωuTrα1+iωuTrα(28) can be used to find a closed-form approximation of the eigenvalue ω~d. Of course, -ω~d∗ 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 massme, viscous damping ce, and stiffness ke, are required for its definition. As is well known, the FRF of that model will have the following form:(34)Heω=1-meω2+ceiω+ke≡1me1-ω2+2iωωeζe+ωe2≡1me1Leωwhere theequivalent undamped natural frequencyωe=ke/me and itsequivalent damping ratioζe=ce/2meωe have been introduced. Our first condition is to impose the fact that the poles of He(ω) are the same as those of the original FRF, H(ω). Consequently, the main feature of the approached function Le(ω) is to take the form of a second-order polynomial with roots ω~d and -ω~d∗. Thus,(35)Leω=-ω-ω~dω+ω~d∗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:(36)ω~d=a+biωu,ω~d∗=a-biωuwhere a,b>0 are two nondimensional parameters. It is worth highlighting that ω~d=a+ibωu is analytically available from (28) and therefore both the real part a=Rω~d and the imaginary part Iω~d are closed-form expressions explicitly defined as functions of the parameters of the original model. Mathematically, they can be expressed as certain functions fa and fb such that(37)a=12ω~d+ω~d∗ωu=fam,k,α,c,Trb=12iω~d-ω~d∗ωu=fbm,k,α,c,TrIntroducing (36) into (35) and expanding the product result in(38)Leω=-ω2+2iωbωu+a2+b2ωu2=-ω2+2iωba2+b2a2+b2ωu+a2+b2ωu2=-ω2+2iωζeωe+ωe2Identifying the last two expressions, the damping factor and the natural frequency of the equivalent viscous model can be obtained as functions of a,b and ωu.(39)ζe=ba2+b2,ωe=a2+b2ωuSecondly, 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 ζe and natural frequency ωe, the peak of the absolute value of the FRF is located in the abscissa ωmax=ωe1-2ζe2. Due to the relevance of the range around this point in the response, it is reasonable to calculate the mass me under the condition that both FRFs have the same value at ω=ωmax. Mathematically, this can be expressed by the equality between both models.(40)Hωmax=1m1Lωmax=1me1Leωmax=HeωmaxTherefore, the mass of the equivalent model can be expressed by(41)me=LωmaxLeωmaxmUsing (39), ωmax results:(42)ωmax=a2-b2ωu=a1-ba2ωu≡a¯ωuSubstituting this expression in Lωmax and Leωmax and after some straight operations, a more compact expression for me can be deduced:(43)me=φa¯ωu-a¯22abmTo sum up, it has been proven that the three equivalent parameters, me,ce,ke, are directly related to the parameters of the original model m,k,α, and Tr through the relations(44)me=λ2abmce=2meζeωe=λamωuke=meωe2=λ2ab+bakand using (37). Above, the notation λ=φ(a¯ωu)-a¯2 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 matrixHω 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 ω.(45)v^ω=-ω2M+φωK-1F^ω≡HωF^ω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 K and M. Hence, a matrix Ψ=diagψ1,…,ψn∈Rn×n exists, such that(46)ΨTMΨ=M′,ΨTKΨ=K′where M′=diag[mj′] and K′=diag[kj′] are two diagonal matrices whose entries are related to the natural frequencies Λ=diag[ωj2] through the equation K′=ΛM′. The mass-normalization of each eigenvector ψj shows as a result the matrix Φ=Ψ(M′)-1/2 previously presented in (19). Thus, the receptance matrix is obtained from the following matrix product:(47)Hω=Ψ·diag1-mj′ω2+φωkj′·ΨT≡ΨTωΨTEach entry of the main diagonal in Tω is the FRF of the jth modal coordinate, Tjω. The resulting FRF has the same form as that of a SDOF system defined in (32) with mass mj′ and stiffness kj′. Thus, introducing again the function Ljω, which only depends on the parameters of the damping model α,c,Tr and on the jth undamped frequency, the jth modal FRF leads to(48)Tjω=1mj′1-ω2+φωωj2≡1mj′1LjωThe set of complex frequencies ω~j=aj+bjiωjj=1n are computed using the above approach (28). Thus, as a simple extension of the method described for the SDOF systems, the function Tjω can be approximated by(49)Tjω≈Tejω=1-mej′ω2+cej′iω+kej′=1mej′1-ω2+2iωωejζej+ωej2where the equivalent parameters mej′, cej′, and kej′ are given by(50)mej′=λj2ajbjmj′,cej′=λjajmj′ωj,kej′=12ajbj+bjajλjkj′(51)ζej=bjaj2+bj2,ωej=aj2+bj2ωj,1≤j≤nwith λj=φ(a¯jωj)-a¯j2 and a¯j=aj1-bj2/aj2. These equations can be put in matrix form by(52)Me′=M′·Π′,Ke′=K′·Γ′,Ce′=2Me′ΛeZewhere Λe=diagωej2 and Ze=diagζej and the two diagonal matrices, Π′ and Γ′, have been introduced, which are working as transformation matrices in the modal space:(53)Π′=diagλj2ajbj,Γ′=diag12ajbj+bjajλj.Therefore, the approximation of Tω is given by(54)Tω≈Teω=-ω2Me′+iωCe′+Ke′-1Furthermore, the receptance matrix Hω of the MDOF system can be estimated with(55)Hω≈Heω=ΨTeωΨTThe above expression is just the classical form of a receptance matrix corresponding to a linear-viscous model. Now, three matricesMe, Ce, and Ke can be defined so that they verify the following conditions:(56)Me′=ΨTMeΨCe′=ΨTCeΨKe′=ΨTKeΨFinally, the proposed viscous approach for MDOF viscoelastic systems is defined as the linear system (57)Mev¨+Cev˙+Kev=Ftand has He(ω) 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 Me and Ke. Thus, starting from their definition and using (46) and (52), we have(58)Me=Ψ-TM′Π′Ψ-1=Ψ-TΨTMΨΠ′Ψ-1≡MΠKe=Ψ-T·K′Γ′·Ψ-1=Ψ-TΨT·K·ΨΠ′Ψ-1≡K·ΓCe=2MΠΨΛeZeΨ-1where Ψ-T=ΨT-1=Ψ-1T denotes the inverse transpose of Ψ. Above, the following matrices have been introduced:(59)Π=ΨΠ′Ψ-1,Γ=ΨΓ′Ψ-1which can be interpreted as transformation matrices affecting the original matrices M and K to provide as an output their equivalent matrices Me and Ke, 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, He(ω), is a good approximation of the original model one, H(ω). 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ηm given in (10).


## 4.1. Complex Eigenvalues

Let us consider the eigenvalue problem corresponding to the homogeneous system in (15).(18)-ω2M+φωKv^=0^In general, the number of eigenvalues of this problem can be greater than 2n, say 2n+p, depending on the nature of the viscoelastic model [45]. The p additional roots are namednonviscous modes and are in general overcritically damped. In this paper, we only consider the 2nelastic 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 ω~j,-ω~j∗j=1n, 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ϕj∈Rn the jth mass-normalized jth eigenvector and ωj is its undamped natural frequency. Arranging eigensolutions in the matrices Φ=[ϕ1,…,ϕn] and Λ=diag[ω12,…,ωn2], both within Rn×n, the following relations are well known from the classical modal analysis:(19)ΦTMΦ=In,ΦTKΦ=Λwhere In is the nth-order identity matrix. The change of variable v^=Φy^ transforms in a diagonal system the eigenvalue problem of (18). Hence, the complex eigenvalues are the solution of the equation(20)-ω2In+φωΛ=∏j=1n-ω2+φωωj2≡∏j=1nLjω=0In particular, the jth complex eigenfrequency must verify the equation(21)Ljω=-ω2+φωωj2=0,j=1,…,nwhere ωj is the jth 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 (s=iω):(22)djs=s2+sgs+ωj2=0The key idea is the expansion of dj(s) around s0=-ξωj±ωj1-ξ2, where ξ=lims→0G(s)/2. It can be easily shown that our equation (21) can be expressed as the form (22), taking a viscoelastic function defined by(23)gs=φs-1swhere φs is the function defined by (7) but changing iω→s. However, in this particular case, a singularity appears when the previous limit is evaluated. Actually, applying L'Hôpital's rule results in(24)lims→0⁡gs=lims→0⁡φs-1s=lims→0⁡φ′s=lims→0⁡1sTr1-ααc-11+sTrα=+∞assuming that 0<α<1 and c>1. Therefore, it is not suitable to apply the method just as it is described in [34]. However, the expansion of Ljω is always possible provided that another valid initial point can be found. In this work, the following alternative initial point is proposed:(25)s0=iωjφωjThis approximation comes from solving the characteristic equation after evaluating the viscoelastic function at the jth undamped natural frequency; that is,(26)-ω2+φωjωj2=0⟶ω~j0=+ωjφωj∈CThe 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 ω~j0 are expected as well. As a consequence of that, expanding Lj(ω~j) in its Taylor series up to the second order should be a good approximation for the complex frequencies. Thus, if ω~j is a root, there exists a number ε∈C, such that ω~j=ω~j0+ε. Thus, a second-order polynomial results:(27)Ljω~j0+ε≈Ljω~j0+ε∂Lj∂ωω=ω~j0+12ε2∂2Lj∂ω2ω=ω~j0=0with the following complex roots:(28)ω~j≈ω~j0+2ω~j0-ωj2φ0′±ωj2φ0′-2ω~j02-2ωj2φ0-ω~j0ωj2φ0″-2φ0″ωj2-2where(29)φ0=φω~j0,φ0′=∂φ∂ωω=ω~j0,φ0′′=∂2φ∂ω2ω=ω~j0Among both roots of (28), that one with positive real part is chosen, assuming that such value will be closer to ω~j0. 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). Figure1 shows schematically a SDOF system with a mass m attached to a fixed point through a viscoelastic constraint modeled with fractional derivatives.Figure 1
The original single-degree-of-freedom model with viscoelastic behavior based on the fractional derivatives and its equivalent linear viscous model.Dα(•) denotes the α-order fractional derivative.LetR(t) and F(t) denote the force reaction and the time-dependent external force, respectively. The relationship between R(t) and the degree of freedom u(t) can be expressed as (see (4))(30)R+TrαdαRdtα=ku+cTrαdαudtαwhere α, c, and Tr are the parameters of the viscoelastic model based on fractional derivatives. Applying the Fourier transform to this expression, it follows that R^(ω)=k^(ω)u^(ω), where k^(ω)=φ(ω)k denotes the complex stiffness [2]. Assuming an excitation force F(t), the motion equation in frequency domain can be expressed as(31)-mω2+φωku^=F^ωand we introduce H(ω) as the frequency response function (FRF) of the system:(32)u^=HωF^ω,Hω=1m1-ω2+φωωu2≡1m1Lωwhere ωu=k/m is the undamped natural frequency of the system. Now, using as an initial point(33)ω~d0=ωuφωu=ωu1+ciωuTrα1+iωuTrα(28) can be used to find a closed-form approximation of the eigenvalue ω~d. Of course, -ω~d∗ 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 massme, viscous damping ce, and stiffness ke, are required for its definition. As is well known, the FRF of that model will have the following form:(34)Heω=1-meω2+ceiω+ke≡1me1-ω2+2iωωeζe+ωe2≡1me1Leωwhere theequivalent undamped natural frequencyωe=ke/me and itsequivalent damping ratioζe=ce/2meωe have been introduced. Our first condition is to impose the fact that the poles of He(ω) are the same as those of the original FRF, H(ω). Consequently, the main feature of the approached function Le(ω) is to take the form of a second-order polynomial with roots ω~d and -ω~d∗. Thus,(35)Leω=-ω-ω~dω+ω~d∗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:(36)ω~d=a+biωu,ω~d∗=a-biωuwhere a,b>0 are two nondimensional parameters. It is worth highlighting that ω~d=a+ibωu is analytically available from (28) and therefore both the real part a=Rω~d and the imaginary part Iω~d are closed-form expressions explicitly defined as functions of the parameters of the original model. Mathematically, they can be expressed as certain functions fa and fb such that(37)a=12ω~d+ω~d∗ωu=fam,k,α,c,Trb=12iω~d-ω~d∗ωu=fbm,k,α,c,TrIntroducing (36) into (35) and expanding the product result in(38)Leω=-ω2+2iωbωu+a2+b2ωu2=-ω2+2iωba2+b2a2+b2ωu+a2+b2ωu2=-ω2+2iωζeωe+ωe2Identifying the last two expressions, the damping factor and the natural frequency of the equivalent viscous model can be obtained as functions of a,b and ωu.(39)ζe=ba2+b2,ωe=a2+b2ωuSecondly, 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 ζe and natural frequency ωe, the peak of the absolute value of the FRF is located in the abscissa ωmax=ωe1-2ζe2. Due to the relevance of the range around this point in the response, it is reasonable to calculate the mass me under the condition that both FRFs have the same value at ω=ωmax. Mathematically, this can be expressed by the equality between both models.(40)Hωmax=1m1Lωmax=1me1Leωmax=HeωmaxTherefore, the mass of the equivalent model can be expressed by(41)me=LωmaxLeωmaxmUsing (39), ωmax results:(42)ωmax=a2-b2ωu=a1-ba2ωu≡a¯ωuSubstituting this expression in Lωmax and Leωmax and after some straight operations, a more compact expression for me can be deduced:(43)me=φa¯ωu-a¯22abmTo sum up, it has been proven that the three equivalent parameters, me,ce,ke, are directly related to the parameters of the original model m,k,α, and Tr through the relations(44)me=λ2abmce=2meζeωe=λamωuke=meωe2=λ2ab+bakand using (37). Above, the notation λ=φ(a¯ωu)-a¯2 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 matrixHω 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 ω.(45)v^ω=-ω2M+φωK-1F^ω≡HωF^ω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 K and M. Hence, a matrix Ψ=diagψ1,…,ψn∈Rn×n exists, such that(46)ΨTMΨ=M′,ΨTKΨ=K′where M′=diag[mj′] and K′=diag[kj′] are two diagonal matrices whose entries are related to the natural frequencies Λ=diag[ωj2] through the equation K′=ΛM′. The mass-normalization of each eigenvector ψj shows as a result the matrix Φ=Ψ(M′)-1/2 previously presented in (19). Thus, the receptance matrix is obtained from the following matrix product:(47)Hω=Ψ·diag1-mj′ω2+φωkj′·ΨT≡ΨTωΨTEach entry of the main diagonal in Tω is the FRF of the jth modal coordinate, Tjω. The resulting FRF has the same form as that of a SDOF system defined in (32) with mass mj′ and stiffness kj′. Thus, introducing again the function Ljω, which only depends on the parameters of the damping model α,c,Tr and on the jth undamped frequency, the jth modal FRF leads to(48)Tjω=1mj′1-ω2+φωωj2≡1mj′1LjωThe set of complex frequencies ω~j=aj+bjiωjj=1n are computed using the above approach (28). Thus, as a simple extension of the method described for the SDOF systems, the function Tjω can be approximated by(49)Tjω≈Tejω=1-mej′ω2+cej′iω+kej′=1mej′1-ω2+2iωωejζej+ωej2where the equivalent parameters mej′, cej′, and kej′ are given by(50)mej′=λj2ajbjmj′,cej′=λjajmj′ωj,kej′=12ajbj+bjajλjkj′(51)ζej=bjaj2+bj2,ωej=aj2+bj2ωj,1≤j≤nwith λj=φ(a¯jωj)-a¯j2 and a¯j=aj1-bj2/aj2. These equations can be put in matrix form by(52)Me′=M′·Π′,Ke′=K′·Γ′,Ce′=2Me′ΛeZewhere Λe=diagωej2 and Ze=diagζej and the two diagonal matrices, Π′ and Γ′, have been introduced, which are working as transformation matrices in the modal space:(53)Π′=diagλj2ajbj,Γ′=diag12ajbj+bjajλj.Therefore, the approximation of Tω is given by(54)Tω≈Teω=-ω2Me′+iωCe′+Ke′-1Furthermore, the receptance matrix Hω of the MDOF system can be estimated with(55)Hω≈Heω=ΨTeωΨTThe above expression is just the classical form of a receptance matrix corresponding to a linear-viscous model. Now, three matricesMe, Ce, and Ke can be defined so that they verify the following conditions:(56)Me′=ΨTMeΨCe′=ΨTCeΨKe′=ΨTKeΨFinally, the proposed viscous approach for MDOF viscoelastic systems is defined as the linear system (57)Mev¨+Cev˙+Kev=Ftand has He(ω) 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 Me and Ke. Thus, starting from their definition and using (46) and (52), we have(58)Me=Ψ-TM′Π′Ψ-1=Ψ-TΨTMΨΠ′Ψ-1≡MΠKe=Ψ-T·K′Γ′·Ψ-1=Ψ-TΨT·K·ΨΠ′Ψ-1≡K·ΓCe=2MΠΨΛeZeΨ-1where Ψ-T=ΨT-1=Ψ-1T denotes the inverse transpose of Ψ. Above, the following matrices have been introduced:(59)Π=ΨΠ′Ψ-1,Γ=ΨΓ′Ψ-1which can be interpreted as transformation matrices affecting the original matrices M and K to provide as an output their equivalent matrices Me and Ke, 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, He(ω), is a good approximation of the original model one, H(ω). 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ηm 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. Figure1 shows schematically a SDOF model with a mass m=1 kg attached to a fixed point through a linear spring with stiffness k=5 kN/m. The viscoelastic model directly depends on the parameters α, c, and Tr and, indirectly, on the loss factor peak ηm (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).Table 1
Example 1: single-degree-of-freedom system. Parameters of the proposed equivalent viscously damped system. Mass and stiffness of the SDOF model:m=1 kg and k=5 kN/m.
Viscoelastic properties Equivalent viscous parameters α c T r ( × 1 0 - 3 s ) η m ζ e ω e / ω u m e / m k e / k CASE 1 0,80 1,014 14,14 0,005 0,0025 1,004 0,998 1,007 CASE 2 0,70 1,165 2,83 0,050 0,0168 1,015 0,989 1,030 CASE 3 0,60 1,981 0,71 0,250 0,0513 1,052 0,971 1,106 CASE 4 0,50 5,850 0,28 1,000 0,1518 1,256 0,911 1,578The results presented in Table1 show that the higher the viscoelasticity of the system is, the larger the differences between the equivalent parameters and the real ones, k and m, 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, ηm.Figures2(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 ηm 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, ζe and ωe, 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 me 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, ω=0. Thus(60)H0He0=memLe0L0=meωe2mωu2=kekAccording to Table 1, the value ke/k increases with ηm. 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.Figure 2
Example 1: single-degree-of-freedom system. Frequency response functions of viscoelastic (exact) and viscous (proposed) system. The four cases covered in Table1 are represented from plots (a) to (d), respectively.
(a)
CASE 1,ηm=0.005 (b)
CASE 2,ηm=0.050 (c)
CASE 3,ηm=0.250 (d)
CASE 4,ηm=1.000In 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 as(61)ϵω=Hω-HeωHωhas been plotted for the four cases of study in Figure 3. The curves show a clear correlation between the loss factor peak ηm 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.Figure 3
Example 1: single-degree-of-freedom system. Relative error (%) between frequency response functions of the viscoelastic model (exact) and the viscous approach (proposed) for the four cases shown in Table1.


### 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 Sections2 and 3 is used. The finite element model is formed by ne=24 finite elements. Two degrees of freedom are associated with each node: the vertical displacement w and the rotation θz. 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.Table 2
Example 2: viscoelastic 1D structures. Viscoelastic material and beam model mechanical properties. The reference frequency isω0=EIy/ρAL4, where A and Iy are, respectively, cross-sectional area and moment of inertia.
Density Young’s modulus Shear modulus Length Ref. frequency Section Elements ρ (t/m3) E (GPa) G (GPa) L (m) ω 0 (Hz) b × h (mm) n e 1,25 10 3,7 1,00 13,00 100×50 24Figure 4
Example 2: viscoelastic 1D structures. Location of the excitation force and the response for the computation of the frequency response function (FRF).As for the first example, the viscoelastic parameters will be presented and ordered in Table3, considering the four numerical cases, according to the level of viscoelasticity, measured by the loss factor peak, ηm. As a first step, the set of complex frequencies ω~j=aj+bjiωj,1≤j≤n, need to be calculated. They will be approximated using the proposed approach given by (28). It is expected that both the real part aj and the imaginary part bj 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 eRe and eIm denote the maximum error in real imaginary part, respectively. Mathematically, they can be found by(62)eRe=max1≤j≤n⁡Rω~j,approx-Rω~j,exactRω~j,exact,eIm=max1≤j≤n⁡Iω~j,approx-Iω~j,exactIω~j,exactwhere n 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.Table 3
Example 2: viscoelastic 1D structures. Damping parameters of the viscoelastic model based on fractional derivatives and numerical cases defined for the beam model. Last three columns show the maximum error in the real part,eRe, and imaginary part, eIm, in the computation of the complex frequencies.
Viscoelastic properties ErrorR(ω~j) ErrorI(ω~j) α c T r ( × 1 0 - 3 s ) η m e R e (%) e I m (%) CASE 1 0,80 1,014 11,24 0,005 2,71⋅10-14 2,03⋅10-10 CASE 2 0,70 1,165 1,35 0,050 8,78⋅10-11 1,59⋅10-8 CASE 3 0,60 1,981 0,34 0,250 6,62⋅10-7 6,31⋅10-5 CASE 4 0,50 5,850 0,07 1,000 5,85⋅10-5 5,80⋅10-3Frequency response function (FRF) of the exact model,Hij(ω), and its viscous approximation, He,ij(ω), 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 Hijω as a curve in the complex plane, controlled by the parameter ω. When the curve Hijω 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 Hijω 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).Figure 5
Example 2: viscoelastic 1D structures. Frequency response functions of viscoelastic (exact) and viscous (proposed) system. The four cases covered in Table1 are represented from plots (a) to (d), respectively.
(a)
CASE 1,ηm=0.005 (b)
CASE 2,ηm=0.050 (c)
CASE 3,ηm=0.250 (d)
CASE 4,ηm=1.000Figure 6
Example 2: viscoelastic 1D structures. Relative error (%) between frequency response functions of the viscoelastic model (exact) and the viscous approach (proposed) for the four cases shown in Table3.As it has been seen before, it is interesting to compare the four damping cases listed in Table3 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 (ω=0), arise from the difference between FRFs at ω=0. In fact, at ω=0, it yields H(0)=K≠He(0)=Ke. In the paper, it is also demonstrated that differences between matrices K and Ke 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.


## 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. Figure1 shows schematically a SDOF model with a mass m=1 kg attached to a fixed point through a linear spring with stiffness k=5 kN/m. The viscoelastic model directly depends on the parameters α, c, and Tr and, indirectly, on the loss factor peak ηm (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).Table 1
Example 1: single-degree-of-freedom system. Parameters of the proposed equivalent viscously damped system. Mass and stiffness of the SDOF model:m=1 kg and k=5 kN/m.
Viscoelastic properties Equivalent viscous parameters α c T r ( × 1 0 - 3 s ) η m ζ e ω e / ω u m e / m k e / k CASE 1 0,80 1,014 14,14 0,005 0,0025 1,004 0,998 1,007 CASE 2 0,70 1,165 2,83 0,050 0,0168 1,015 0,989 1,030 CASE 3 0,60 1,981 0,71 0,250 0,0513 1,052 0,971 1,106 CASE 4 0,50 5,850 0,28 1,000 0,1518 1,256 0,911 1,578The results presented in Table1 show that the higher the viscoelasticity of the system is, the larger the differences between the equivalent parameters and the real ones, k and m, 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, ηm.Figures2(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 ηm 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, ζe and ωe, 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 me 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, ω=0. Thus(60)H0He0=memLe0L0=meωe2mωu2=kekAccording to Table 1, the value ke/k increases with ηm. 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.Figure 2
Example 1: single-degree-of-freedom system. Frequency response functions of viscoelastic (exact) and viscous (proposed) system. The four cases covered in Table1 are represented from plots (a) to (d), respectively.
(a)
CASE 1,ηm=0.005 (b)
CASE 2,ηm=0.050 (c)
CASE 3,ηm=0.250 (d)
CASE 4,ηm=1.000In 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 as(61)ϵω=Hω-HeωHωhas been plotted for the four cases of study in Figure 3. The curves show a clear correlation between the loss factor peak ηm 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.Figure 3
Example 1: single-degree-of-freedom system. Relative error (%) between frequency response functions of the viscoelastic model (exact) and the viscous approach (proposed) for the four cases shown in Table1.


## 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 Sections2 and 3 is used. The finite element model is formed by ne=24 finite elements. Two degrees of freedom are associated with each node: the vertical displacement w and the rotation θz. 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.Table 2
Example 2: viscoelastic 1D structures. Viscoelastic material and beam model mechanical properties. The reference frequency isω0=EIy/ρAL4, where A and Iy are, respectively, cross-sectional area and moment of inertia.
Density Young’s modulus Shear modulus Length Ref. frequency Section Elements ρ (t/m3) E (GPa) G (GPa) L (m) ω 0 (Hz) b × h (mm) n e 1,25 10 3,7 1,00 13,00 100×50 24Figure 4
Example 2: viscoelastic 1D structures. Location of the excitation force and the response for the computation of the frequency response function (FRF).As for the first example, the viscoelastic parameters will be presented and ordered in Table3, considering the four numerical cases, according to the level of viscoelasticity, measured by the loss factor peak, ηm. As a first step, the set of complex frequencies ω~j=aj+bjiωj,1≤j≤n, need to be calculated. They will be approximated using the proposed approach given by (28). It is expected that both the real part aj and the imaginary part bj 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 eRe and eIm denote the maximum error in real imaginary part, respectively. Mathematically, they can be found by(62)eRe=max1≤j≤n⁡Rω~j,approx-Rω~j,exactRω~j,exact,eIm=max1≤j≤n⁡Iω~j,approx-Iω~j,exactIω~j,exactwhere n 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.Table 3
Example 2: viscoelastic 1D structures. Damping parameters of the viscoelastic model based on fractional derivatives and numerical cases defined for the beam model. Last three columns show the maximum error in the real part,eRe, and imaginary part, eIm, in the computation of the complex frequencies.
Viscoelastic properties ErrorR(ω~j) ErrorI(ω~j) α c T r ( × 1 0 - 3 s ) η m e R e (%) e I m (%) CASE 1 0,80 1,014 11,24 0,005 2,71⋅10-14 2,03⋅10-10 CASE 2 0,70 1,165 1,35 0,050 8,78⋅10-11 1,59⋅10-8 CASE 3 0,60 1,981 0,34 0,250 6,62⋅10-7 6,31⋅10-5 CASE 4 0,50 5,850 0,07 1,000 5,85⋅10-5 5,80⋅10-3Frequency response function (FRF) of the exact model,Hij(ω), and its viscous approximation, He,ij(ω), 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 Hijω as a curve in the complex plane, controlled by the parameter ω. When the curve Hijω 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 Hijω 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).Figure 5
Example 2: viscoelastic 1D structures. Frequency response functions of viscoelastic (exact) and viscous (proposed) system. The four cases covered in Table1 are represented from plots (a) to (d), respectively.
(a)
CASE 1,ηm=0.005 (b)
CASE 2,ηm=0.050 (c)
CASE 3,ηm=0.250 (d)
CASE 4,ηm=1.000Figure 6
Example 2: viscoelastic 1D structures. Relative error (%) between frequency response functions of the viscoelastic model (exact) and the viscous approach (proposed) for the four cases shown in Table3.As it has been seen before, it is interesting to compare the four damping cases listed in Table3 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 (ω=0), arise from the difference between FRFs at ω=0. In fact, at ω=0, it yields H(0)=K≠He(0)=Ke. In the paper, it is also demonstrated that differences between matrices K and Ke 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.



---

*Source: 5957831-2018-06-21.xml*