Analysis of Nonviscous Oscillators Based on the Damping Model Perturbation

A novel numerical approach to compute the eigenvalues of linear viscoelastic oscillators is developed.The dissipative forces of these systems are characterized by convolution integrals with kernel functions, which in turn contain a set of damping parameters. The free-motion characteristic equation defines implicitly the eigenvalues as functions of such parameters. After choosing one of them as independent variable, the key idea of the current paper is to obtain a differential equation whose solution can be considered, under certain conditions, a good approximation. The method is validated with several numerical examples related to damping models based on exponential kernels, on fractional derivatives, and on the well-known viscous model. Taylor series expansions up to the second order are obtained and in addition analytical solutions for the viscous model are achieved. The numerical results are very close to the exact ones for light and medium levels of damping and also very good for high levels if the chosen parameter is close to initial values that are defined for every case.


Introduction
The importance of mathematical modeling for structural dynamic systems in fields such as Aerospace or Civil Engineering is well known.One of the most studied issues in the last 50 years has been the search for robust and efficient models for energy dissipation mechanisms.In this sense, the linear viscoelastic models have been profusely used, covering a wide range of structural systems.A large number of analysis methods for these models depend on the solution of the associated eigenvalue problem.This dependency has motivated us to investigate efficient numerical methods to compute the complex eigenvalues of linear viscoelastic oscillators.
The most general case of a viscoelastic system with one degree-of-freedom, say (), can be modeled with a single mass, say , attached to a fixed point by a viscoelastic constraint.Figure 1 shows the schematic configuration mass-spring-viscoelastic damper and the corresponding free body diagram with the applied forces of the mass.Hence, the reaction force produced in the mass, say (), is related with the displacement by where G() is the dissipative kernel or damping function, characteristic of the considered viscoelastic model, and  is the constant of the linear-elastic spring.Notice that the dissipative force depends on the past history of velocities via convolution integral over the kernel.
The motion equation can be deduced from the dynamic equilibrium: where  0 and V 0 are the initial position and velocity of the mass and () is the external, time-dependent, applied force.With reference to the mathematical form of G(), several works have been published, proposing different expressions basically in the Laplace domain.The most popular is due to Biot [1] and is called multiexponential model, which assumes that the viscoelastic function is a linear combination of exponential functions.Also, the models based on the fractional derivative proposed by Bagley and Torvik [2,3] have been widely used.Other viscoelastic models of special interest, again based on state-space methods, are the GHM approach of Golla, Hughes, and McTavish [4,5] or the Anelastic Displacement Field of Lesieutre and Mingori [6].Reference [4] gave the required conditions of G() to describe a real dissipative motion.Adhikari and Woodhouse [7,8] proposed viscoelastic functions in the context of damping identification in dynamic systems.
The main objective of the present work is to propose a numerical method to solve the characteristic equation associated with the nonviscously damped oscillator described by (2).As known, the approximation embraces not only physical systems with only one degree-of-freedom but also multiple-degree-of-freedom systems with proportional or lightly nonproportional nature, that is, those in which the dynamic matrices become diagonal in the modal space.The Laplace transform of the free-motion equation leads to a nonlinear eigenvalue problem in the Laplace variable .The form of this problem will depend on the nature of the viscoelastic function transform, named () = L{G()}.Several methods to solve the general nonlinear eigenvalue problem exist in the bibliography.Yang [9] and Singh and Ram [10] proposed some methods based on Taylor series expansion of the transcendental matrices combined with Newton's eigenvalue iteration method.Williams and Kennedy [11] developed a numerical procedure based on the parabolic interpolation of the determinant that governs the eigenvalue problem.Daya and Potier-Ferry [12] used asymptotic numerical techniques to determine the natural frequencies and the loss factors of viscoelastic damped sandwich structures.Voss [13,14] proposed two methods based on the shift-and-invert Arnoldi's technique and on the Jacobi-Davidson method, respectively.In both references, an iterative process is used to compute a few eigenvalues that lie close to a given point.An important subset of viscoelastic models, in which () has a rational expression, have been published.In them, the nonlinear eigenvalue problem can be transformed into a linear one, introducing new internal variables and using state-space approaches; see for instance the works of Muravyov and Menon [15][16][17].However, these methods become computationally expensive for large systems.In addition, the physical insight of the problem can also be lost due to the introduction of the new variables.Recently, Adhikari and Pascual [18,19] have published efficient iterative methods based on the Taylor series expansion of ().Lázaro et al. [20] proposed a recursive scheme based on the fixed-point iteration which always converge to the complex eigenvalues.
Several authors have studied the eigensolutions as functions of the parameters involved in the dynamic matrices.The main results are related to numerical methods for the computation of eigenvalues and eigenvector derivatives with respect to certain design parameter.Fox and Kapoor [21] published an important work in which analytical solutions of the derivatives were given.However, in their results, the complete set of eigenvectors of the modal space are involved, and again the computations become expensive for largeorder systems.Nelson [22] gave an alternative and efficient method to obtain the derivatives requiring only the eigenvector and eigenvalue of interest.Wang [23] and later Zhang and Zerva [24] introduced improved methods in truncated systems.Murthy and Haftka [25] published a survey of methods for the calculation of the derivatives applied to general non-Hermitian matrix problems.The same authors in [26] studied the computational efficiency of different approximations based on eigenvalue derivatives, generalized Rayleigh quotient, and the trace theorem.The generalization of the eigensolution derivatives for nonviscous systems was studied by Adhikari [27,28].Cortés and Elejebarrieta [29,30] used Adhikari's solutions in an iterative numerical method to compute eigensolutions which are applicable even to highly damped systems.
The current paper is aimed at numerically solving the eigenvalues of the nonlinear eigenvalue problem of viscoelastic oscillators.As known, a dynamic model is governed by inertial masses, linear rigidities, and several damping parameters that are included in the viscoelastic function.The newly developed eigenvalues will be functions of a single damping parameter, called valid parameter, leaving the others fixed.The key idea is to transform the modal characteristic equation into a differential one with separated variables.The methodology is applied to widely used viscoelastic systems and also illustrated with numerical examples to validate the theoretical results and to show the limits of application.

Dynamics of Single Degree-of-Freedom Viscoelastic Systems
The characteristic equation of the eigenvalue problem can be obtained trying solutions with the form () =  0   in the free-motion equation obtained from (2).Consider Here, as mentioned, () is the Laplace transform of the kernel function G().Note that the eigenvalue problem of proportional or lightly nonproportional MDOF systems of  degrees of freedom can be reduced to solve just  decoupling characteristic equations; therefore, the method described in this paper can also be applied to these MDOF dynamic systems.Adhikari [31] gave the necessary and sufficient conditions under which nonviscous systems have proportional damping.The previous equation can be mass-normalized to The damping function Γ() = ()/ and the natural frequency of the undamped system   = √/ have been introduced.Golla and Hughes [4] gave the necessary conditions that () has to fulfil in order to define a real dissipative motion.Assuming that () verifies these conditions, it can be assured [32,33] that the characteristic equation ( 4) has 2+ eigenvalues, with  ≥ 0. These can be arranged in a set with the form where  0 ,  * 0 are a pair of complex conjugate roots which induce oscillatory motion.The rest {  }  =1 are named nonviscous eigenvalues [34] and are negative real numbers associated with overcritical damped responses.
According to Adhikari [34], if the set of eigenvalues is available, the complete solution of the integrodifferential equation ( 2) can be expressed in terms of these eigenvalues in the form where the function and coefficient are Due to the relevance of the eigenvalues in the final solution, the availability of efficient numerical computation tools is important.We propose a new method to compute the roots of (3), based on the perturbation of the damping parameters.In the next section, the theoretical foundations of the method are introduced.

Perturbation Method of the Damping Model
3.1.Fundamentals of the Method.As explained in Section 2, the eigenvalues are the roots of the characteristic equation given by (4).The key idea of the method is to consider that the function Γ = Γ(,  1 , . . .,   ) depends not only on the variable  but also on a set of damping parameters { 1 , . . .,   } which controls the viscoelastic behavior of the system.
In the general case, any eigenvalue  = ( 1 , . . .,   ) can be considered as function of the damping parameters.Let us consider one of these parameters as variable, while the rest remain fixed.Without loss of generality, it can be assumed that this is  1 .Denoting as  =  1 , let us assume now that there exists a particular value (or initial point) of , say  0 , so that the characteristic equation can be solved analytically in terms of  2 ,  3 , . . .,   , obtaining the related initial eigenvalue  0 : Under these conditions, the parameter  =  1 is named valid parameter for the purposes of the current paper.Somehow, this parameter will be used to perturb (4) and therefore we need a particular initial value, which has been denoted by  0 .The method is constructed assuming that the eigenvalue at  =  0 , say  0 , can be analytically available (and consequently also its derivatives will be shown later).Therefore, we can take as independent variables those parameters which, evaluated at certain particular value, allow reducing (4) to another one with analytical solution.Along the numerical examples of the present paper, different choices for the damping parameters will be carried out in order to show the flexibility of the proposed method.Since the rest of parameters are fixed, the eigenvalue  = () can be considered as single-variable function of .The conditions of existence and uniqueness of such function are given by the implicit function theorem.As known, the equation defines an implicit function  = () around the point Under this hypothesis, it is also guaranteed that there exist (i) an open interval  =] 0 − ℎ/2,  0 + ℎ/2[ with ℎ > 0 and (ii) an open neighborhood around the point  0 in the complex domain  ⊂ C, such that the function () will be unique and continuously differentiable in  and its graph verifies {(,  ()) :  ∈ } = {(, ) ∈  ×  :  (, ) = 0} .(11) Notice that the damping function has been represented directly as Γ(, ) in the implicit function theorem due to the fact that the parameters  2 , . . .,   are assumed to be fixed.The main challenge of this paper is to construct approximated solutions of the function ().In the next subsection, the methodology to obtain the th-order Taylor series of () will be developed.Moreover, this approximation will be used in a further step to improve the solution.

Taylor Expansion of the Eigenvalues.
Assuming that the necessary conditions for the existence of the function  = () are satisfied, the Taylor series expansion up to the th order can be obtained around the initial point  0 , provided that the derivatives evaluated at this point are available.Consider Shock and Vibration with The computation of  () ( 0 ) can be carried out from the characteristic equation (9), in which consecutive derivatives can be taken with respect to .The first-order derivative is calculated as whence solving   results in Using the same procedure to compute the second derivative, the term   is extracted.After some simplifications, this second derivative results in where and now Finding more terms of the series would be possible if higher order derivatives were calculated.However, for the proposes of the present work, only the explicit expressions of the first two derivatives are necessary.Note that both are well defined since the denominators of ( 15) and (17) are never nil due to (10).
As will later be described, in some cases, (12) suffices when the objective is to calculate the eigenvalue at the proximities of point  0 .However, the availability of the solution in a wider interval of the parameter  is sometimes necessary.For these situations, an improved solution is presented in the next subsection.

Improving the Solution: Eigenvalue Differential Equation.
The challenge now is to use the previously obtained results to improve the estimation accuracy of () in a wider range of .To describe the procedure, the expression of the first-order derivative   () given in (15) is rewritten in the form where the new function Ψ(, ) is defined by From a mathematical point of view, (20) can be considered an ordinary differential equation to be complemented with the initial condition  0 = ( 0 ).Obviously, its solution is the same as that of the characteristic equation (8).The fundamental objective of the current method is to transform this ordinary differential equation into another with separated variables.For that, the th-order solution of the Taylor series obtained in (12), namely,   (), will be used to approximate the value of  inside the function Ψ(, ).So Hence, the new function   () depends only on a single variable.The accuracy of the approximation is related to the quality of the assumption made in (22), which in turn is conditioned by the level of viscoelasticity of the system.In other words, it depends on the variability of the viscoelastic function Γ(, ) with respect to ; a viscoelastic system with low viscoelasticity does not present large variations in domain.More details about the study of the viscoelasticity in nonviscous dynamic systems can be found in the work of Adhikari and Woodhouse [35].Here it will be shown that, under certain conditions, the approximation Ψ(, ) ≈ Ψ(  (), ) will produce results close to the exact ones.
With regard to the order of the approximation, in the next subsection, it will be demonstrated that the proposed method improves the error order with respect to that of the Taylor approximation.However, as will be shown with the numerical examples, the order does not ensure a better mean approximation in the wider range of the parameter.Under the aforementioned assumptions, the differential equation can be expressed as whose solution will be denoted by λ ().The subindex indicates the -order Taylor series approximation.An explicit expression can be represented by where the exponential function is defined by When the function λ () is obtained from the first-order approximation of the Taylor series expansion  = 1, the solution will be named Improved Linear Solution (ILS).If the second-order approximation  = 2 is used, the name will be Improved Quadratic Solution (IQS).The functions Ψ(, ) and   () will be called exact and approximated integrand functions, respectively.The next subsection is aimed at the analysis of the truncation error.

Analysis of the Error
Section 3 was focused on the methodology exposition to compute the function λ () as approximation of the exact eigenvalues ().It has already been mentioned that, intuitively, the proposed function will be a good estimation provided that the error between the functions Ψ((), ) and   () is not important.In this subsection more quantitative information is provided, finding a bound of the error | λ () − ()|.The necessary conditions to calculate this bound are directly related with the properties of the viscoelastic function (, ) = Γ(, ) and will be given by three hypotheses: (H2) Let  and  be the neighborhoods of  0 and  0 , respectively; the existence of these sets is assured by the implicit function theorem.Then (, ) is a continuous and differentiable function with respect to the variables ,  up to order  + 1 in the set  × .
In addition, (, ) and its partial derivatives with respect to ,  up to order  + 1 are bounded in the set .
(H3) The function Ψ(, ) defined by ( 21) is Lipschitz continuous with respect to the first variable  in the set , with constant  > 0: The hypothesis (H1) allows assuring the existence of the function () :  →  in the interval  =] 0 − ℎ/2,  0 + ℎ/2[.One of the theses of the implicit function theorem is the existence of the set  ⊂ C in which the image of () is defined.Under this consideration, hypotheses (H2) and (H3) are well established because they are based on the existence of such sets.The main conclusions of the current subsection will be presented in Theorem 1. Previously, two properties are formulated since they will be necessary later in the demonstration of the main theorem.Notice that with the previous nomenclature  (0) () = ().The proof of Properties 1 and 2 is found in Appendices A and B. Theorem 1.Let () be the eigenvalue and let λ () be the approximation calculated from (22) to (25) with  ≥ 1.Assuming that hypotheses (H1), (H2), and (H3) are satisfied, there exist two positive real numbers,   > 0 and  > 0, such that Proof.From ( 20), the exact eigenvalue () verifies the relation where For demonstration purposes, it is necessary to bound the difference norm |  ()−  ()| for any natural number  ≥ 1.
First, using (H3) and the subtraction λ ()−() in Taylor series, the previous norm for  = 1 is where  ∈  defines the truncation residual of the Taylor series up to the order .For the last inequality Property 2 has been used.In general, a bound of the functions () and () can be calculated using Property 1. From the definitions, one directly obtains Second, using ( 30) and ( 31) Now the error can be calculated; using the definitions of both functions, The theorem states that the approximation order of the improved solution is O(ℎ +2 ), while approximating with Taylor series up to the th term gives O(ℎ +1 ).This result allows us to assure that an improvement of the solution in a neighborhood of the initial point exists.
To illustrate the efficiency of the proposed method, the next section describes its implementation for systems with Biot's model-based viscoelastic function.Cases with one exponential kernel and with more than one exponential kernel will be also developed.In addition, the viscously damping systems are analyzed, since the current methodology allows finding analytical solutions.

The Viscous Model.
The viscous model is a special case of viscoelastic behavior whose kernel is proportional to the Dirac delta, G() =  V (), where  V is the viscous damping coefficient.The Laplace transform of the kernel function becomes the -independent function () ≡  V .In general, instead of  V the damping ratio  is commonly used.Both are related by the well-known expression  V = 2  , and consequently the mathematical model depends only on the single parameter  ≡ .The mass-normalized characteristic equation (3) can be rewritten as where now Γ(, ) = 2  .The classical exact solution of (34) is the pair of complex numbers Although the exact solution is here available, the application of the proposed method is of interest due to the availability of analytical solutions for (23) considering Taylor series up to the linear order  1 () or up to the quadratic one  2 ().
Taking  0 = 0 as the initial point, the solution of the characteristic equation corresponds to that of the undamped system, (0) =  0 = ±  .Using now the eigenvalue with positive imaginary part, the derivatives of () can be deduced from ( 15) and ( 17): The Taylor series expansion of the function () up to the second order can be obtained as It can be verified that the three first terms of the Taylor expansion (see (35)) coincide with (37).In view of the theoretical results of the previous section, the improved solutions ILS ( = 1) and IQS ( = 2) can be computed solving the following differential equation: Substituting in (38) the first order approximation of the root  1 () from (37), the ordinary differential equation leads to and a closed expression for ILS is obtained: It is also interesting to deduce the improved solution corresponding to the quadratic approximation  2 () from (37); the differential equation results in and the explicit expression of IQS can also be easily integrated, resulting in Damping ratio,  Both improved solutions together with the exact one are plotted versus the damping ratio  in Figure 2.This ratio is represented in logarithmic scale since the great majority of dynamic systems are lightly damped, with range 0.001 ≤  ≤ 0.2.For this numerical example an undamped natural frequency   = 1rad/s has been chosen.The improved solutions are close or very close to the exact one through almost the complete range of , even for high damped systems.The noticeable differences are for  > 0.6 far from the initial value  0 = 0 as mentioned at the end of Section 4.
Moreover, it can be noted that the improved solutions and the exact one lie along the same complex domain circumference when the parameter  varies.This fact can easily be verified calculating the absolute value of ( 40) and (42): which coincides with the absolute value of the exact root for any damping ratio.

Biot's Exponential
Model.This model was introduced by Biot [1] and has been extensively used for the analysis of viscoelastic systems.It is based on the assumption that dissipative forces depend on the history of the velocities of the degrees of freedom via an exponential kernel with form where  is the relaxation parameter, also called nonviscous parameter.Notice that for high values of  the viscoelastic model tends to the viscous model, mathematically expressed by lim  → ∞ G() =  V ().Therefore, the nonviscous behavior will be controlled by , so that when  ≫   the damping model is viscous with coefficient  V [34,35].Otherwise, when  ≪   the behavior is highly nonviscous.The characteristic equation in the Laplace domain can be written as where () =  V /( + ).The coefficient of the limit viscous model can again be expressed in terms of the damping ratio  V = 2  .Under these conditions, the viscoelastic model can be written as function of the two parameters,  and , as follows: Introducing a new parameter  = 1/ with  > 0, the previous characteristic equation can be rewritten as Equations ( 46) and (47) allow defining the eigenvalues  = () as single-variable functions of either  = ,  = , or  = .Table 1 shows information related with the three cases and the eigenvalues evaluated at the initial values ( 0 ,  0 ).In the present subsection the first two cases will be developed.The third one will be considered for a model with  exponential kernels in Section 5.3.
Case 1 ( = ).The function () is implicitly defined by As an initial point, any root of the equation (, 0) =  2 +  2  = 0 can be taken.Without loss of generality we will use the positive imaginary part  0 = .From ( 15) and ( 17) the values of   (0),   (0) can be computed, obtaining after some operations the expressions Hence, the second-order Taylor polynomial is ILS and IQS, constructed from linear and quadratic Taylor approximations, respectively, can be calculated with (24).For that, the expression of the approximated integrand function must be considered: where   () is defined by (12).As described before, the quality of the approximated eigenvalue function λ () depends on the accuracy of the assumption over the integrand function  5 Ψ((), ) ≈   () given by (22).The method can be implemented calculating the integral given in (25) through numerical quadrature for each parameter value.The results have been represented in Figure 3, where the functions λ1 () and λ2 () have been plotted versus  for an oscillator with natural frequency   = 1 rad/s.Three damping ratios  = {0.01,0.05, 0.10} have been considered for comparison.In addition, the relative error of the exact integrand function Ψ(, ),

Shock and Vibration
is shown.Notice that the exact eigenvalue () from ( 46) is required for the error computation.
Figure 3 shows that, as expected, close to the initial point  = 0, the proposed method eigenvalues are very accurate for all damping levels.Moreover, it can be seen that  2 () <  1 () in the range 0 <  < 1, inequality that, as predicted from the theory of Section 4, implies that IQS is more precise than ILS.On the contrary, when  > 1, the errors are  1 () <  2 () and, in addition, while  2 () increases monotonically,  1 () remains almost constant.Therefore in this range ILS is more accurate and stable than IQS so that the eigenvalues calculated with ILS are in average closer to the exact solution over the considered range of .For a fixed , the higher the damping the poorer the approximation, especially for IQS that completely fails for high values of  and .The divergence of IQS in this situation can be avoided with a different choice of the parameter as will be explained in the next case.
This example suggests that the improved solution with quadratic approximation is not always the best option when the results are required in a wide range.
Case 2 ( =  = 1/).The solution of this case is somewhat equivalent to the expansion of the eigenvalues in their Taylor series around  → ∞ or  = 0.As will be seen, the use of this parameter results in a significant improvement of the results from those of the previous case.The function  = () is now implicitly defined by the equation The complex roots of this equation for  0 = 0 have already been presented in Table 1.As before, (15) and ( 17) are used to analytically calculate the derivatives evaluated at the initial point: With these expressions, the second-order Taylor series expansion of () is defined in the range ,  ≥ 0. Now, ILS and IQS can be calculated with the methodology described in Section 3.3.First, in this particular case, Second, for each value of the parameter , the eigenvalue λ () =  0  −() can be estimated by numerical quadrature of the integral: Figure 4 shows the error and the complex eigenvalues versus  from the exact and the two improved solutions.To compare the results with those of the previous example, the same -axis has been used, taking the inverse of .It is clear that the choice of  produces better results for all  and for the complete range of .The distributions of  1 () and  2 () reverse their relative position but now both tend to zero for high values of .Except close to the origin  = 0 ( → ∞), the estimated eigenvalues are very close to the exact ones.The ILS is again the best solution, although the IQS presents now an almost perfect behavior, even for high damping.

Biot's Multiexponential
Model.This model arises as a natural generalization of Biot's single-exponential model.Using the same nomenclature as that in the previous subsection (see (44)), the viscoelastic kernel function can be written as where   are the relaxation parameters and  V = 2   is, as before, the damping coefficient of the viscous model when   → ∞.The Laplace transform of (58) can be expressed as  (, ,   Directly, the mass-normalized characteristic equation takes the form  negative.The parameters   cannot be used with the proposed method because analytical solutions at the initial point for  kernels are not available.However, if the damping ratio  is used as parameter, the choice of initial value  0 = 0 permits obtaining closed solutions for both the complex eigenvalues and the nonviscous ones. Complex Eigenvalues.The value  0 = 0 produces the eigenvalues of the undamped system in (60); that is,  0 = ±  .Expanding the solution in its Taylor series around this point, the proposed method allows obtaining  = ().The first and second derivatives,   (0) and   (0), can be calculated from ( 15) and (17).After some operations, The second-order Taylor series become function of :  2 () =  0 +   (0) +   (0) 2 /2 ≡  1 () +   (0) 2 /2.As in the previous cases, this approximation can be improved solving the differential equation ( 23), for which The method is implemented for a numerical example with  = 5 kernels.The relaxation parameters are taken as   = {2, 4, 6, 10, 15} rad/s.Three different natural frequencies   = {0.2, 2, 20} rad/s have been considered to validate the results when the system is near-viscous   ≫   or, otherwise, nonviscous   ≈   .
Figure 5 shows the errors of the integrand function and the approximated eigenvalues λ1 () and λ2 () as before.The damping ratio  has again been represented log-scaled.As in the previous case  2 () <  1 (), (several orders of magnitude) ∀ ∈ [0, 1], obtaining better approximations for IQS than those for ILS.For both viscous and nonviscous systems the agreement is very good, except for  ≈ 1 and   = 0.20 rad/s (viscous behavior), where the slope of the imaginary part becomes almost vertical.
Real Eigenvalues.It is expected that the real eigenvalues of (59) will be relatively close to the parameters   in lightly damped systems, as described by Adhikari and Pascual [19].Naming  = () to the nonviscous eigenvalue associated with a particular parameter   , then  0 = lim  → 0 () = −  .Indeed, multiplying (60) by  +   results in Taking limits  → 0 in the previous expression, Since  0 ∈ R for the nonviscous eigenvalues,  2 0 +  2  > 0 and consequently  0 +   = 0. Therefore, the nonviscous eigenvalue associated with   can be expanded around the point  0 = (0) = −  .Without loss of generality, it can always be assumed that () is associated with the first parameter  1 .Thus, for any value  > 0, the Taylor series expansion is now  2 () = − 1 +   (0) +   (0) 2 /2.Since the viscoelastic function and its derivatives are not defined in  = −  (see (59)), the derivatives from ( 15) and ( 17) cannot be calculated.To avoid this problem, (63) with  = 1 will be used.Thus,   (0) and   (0) can be calculated taking derivatives with respect to  in the expression: Evaluating ( 65) at  = 0, and after some operations, the following results are obtained: The computation of the improved solutions from the differential equation (20) requires that () is well defined at the initial point.Since  0 = −  is a pole of () and of ()/, the Taylor series up to the first order and second order will be used to estimate the eigenvalues.
Figure 6 shows results for a numerical example with   = 3rad/s and Biot's model with  = 3 kernels and   = {1, 5, 10} rad/s.As expected, the eigenvalues show very good accuracy for damping range 0 ≤  ≤ 0.20.Better results could be achieved through expansion of the Taylor series up to higher orders.For that, successive derivatives,   (0),  IV (0), . .., from (65) can be extracted and new terms can be added.

Multiple-Degree-of-Freedom Systems.
In this last point, the application of the proposed method for multiple degreesof-freedom (mdof) nonviscously damped systems is presented.First, the derivation of the general expressions will be addressed without formal restrictions on the damping model.Later, a numerical example of a mdof discrete system with a viscoelastic model based on the fractional derivative will be developed and discussed.
Let us consider  degrees-of-freedom (dof) vibrating system.The dofs' time-domain response is represented by a vector u() ∈ R  .With help of the Finite Element Method, the mass and the stiffness matrices of the system, denoted, respectively, by M, K ∈ R × , can be assembled.In general, M, K are symmetric and positive definite and semidefinite, respectively.The viscoelastic damping is introduced in the system assuming that the dissipative forces f  () ∈ R  are proportional to the history of the dofs' velocities via kernel hereditary functions.These functions the entries of  a matrix denoted by G() ∈ R × , also assumed symmetric.Consider Assuming small displacements, the motion equations linearly relate the dynamic balance among inertia M ü, elastic Ku, damping f  (), and external forces f  ().The general form of the system of integrodiferential equations can be written by Under these conditions, the modes of the system can be obtained as the nontrivial solutions of the free-motion problem obtained considering f  () = k 0 = u 0 = 0 in (68).Thus, checking functions of the form u() = u  , we obtain where G() = L{G()} ∈ C × represents the Laplace transform of the viscoelastic kernels and D() : C → C × is the dynamic stiffness matrix.Equation ( 69) is the the nonlinear eigenvalue problem of a viscoelastically damped vibration system.The eigenvalues are then the roots of the equation The complete set of eigenmodes of this problem present 2+ distinct values distributed as where   ,  *  ∈ C, 1 ≤  ≤ , are  conjugate complex pairs, corresponding to the modes with oscillatory nature.The rest,   ∈ R − , 1 ≤  ≤ , are negative real numbers that represent overcritically damped modes called nonviscous since they are a feature of nonviscous systems particularly of those governed by Biot's exponential kernels [36,37].Associated with each eigenvalue there exists an eigenvector: we denote u  , u *  ∈ C  , 1 ≤  ≤  to the complex eigenvectors associated with   ,  *  , and a  ∈ R  , 1 ≤  ≤  to the eigenvector associated with the real eigenvalue   .
In the analysis and search of solutions of the eigenproblem, the proportionality of -dof system, that is, the modal decoupling capability, becomes of special importance.As known, the undamped problem obtained from G() ≡ 0 can always be diagonalized in the modal space of matrices M and K. Denoting by   ∈ R  to the th mass-homogenized real mode (in column) and by   to its associated natural frequency, the classical orthogonal relations can be written as where   is the Kronecker delta.The modal matrix Φ groups in columns  undamped modes and is used to change the dof 's to the modal coordinates by u = Φq.Hence, (69) can be expressed as where is the diagonal squared natural frequency matrix and Γ() = Φ  G()Φ is the damping matrix in the modal space.In general, the latter is not diagonal, although under certain conditions it could become so.The systems with this property are called proportional, since the necessary and sufficient conditions for the modal decoupling are directly related with proportional relationships between the damping matrix G() and the dynamic matrices M, K [31].In some cases of nonproportionality, the matrix Γ(  ) is diagonally dominant but not purely diagonal, fact which is equivalent to assume as true: The system is then said to be lightly nonproportional.This property is commonly assumed in many problems related to nonviscous damping [18,19,34,38] and allows approximating the determinant as product of the terms of its main diagonal as Figure 7: Multiple degrees-of-freedom example: lumped-mass dynamical system with viscoelastic links based on the fractional derivatives.
Hence, the set of eigenvalues can be obtained from the following  decoupled equations: where Γ  () =    G()  .Equation (76) represents the th modal equation and presents the same structure as that of (4) derived for single degree-of-freedom systems.Thus, as done at the first part of this paper, we can assume that some parameter of the damping model is considered variable, say .Consequently, we introduce the notation Γ  (, ) =     G(, )  to refer to the th modal damping function.For certain particular value of this parameter, say  =  0 , (76) admits a closed-form analytic solution, denoted by   ( 0 ).This solution will be taken as our initial point from which we perturb the equation.Therefore, the first-and second-order Taylor series expansions around  =  0 are where the subscripts (⋅) 1 and (⋅) 2 denote the first-and second-order approximations of the th eigenvalue, respectively.Expressions for    and    as function of the damping function can be obtained from (15) and (17).Application of the method includes the definition of the following function for mdof systems and for the mode : as a generalization of that of (21).It has been proved that the Improved Linear and Quadratic Solutions (ILS and IQS, see (24)) defined as present, respectively, one approximation order higher than those of (77) and (78) (see demonstration in Section 4).The computation of the eigenvalue is thus reduced to the calculation of an integral, which can be performed numerically using any quadrature methods available in the literature [39].
To complete the validation of the proposed method on mdof systems, a five-dof discrete system with five masses and viscoelastic links is studied, as shown in Figure 7.The nonviscous constrains (represented in the figure by a linear spring and by a viscoelastic damper) relate reactions and displacements through a four-parameter viscoelastic model based on the fractional derivatives [40].Indexes  = 0, 6 of fixed boundaries are related with zero dof 's so that  0 = u 0 =  6 = u 6 = 0. Thus, the constitutive equations can be written as where Δ  =   −  −1 and , , and   are the parameters of the damping model based on the fractional derivatives, also called storage coefficient, fractional exponent, and relaxation time, respectively.For real materials  > 1, 0 <  < 1,   > 0. The coefficient , is the linear, static rigidity.The kernel function G() is difficult to obtain explicitly and it becomes necessary to appeal to infinite series based functions [41].However, the damping function in the Laplace domain () can easily be calculated by simply applying to the fractional derivatives of (82) the Laplace transform and using its properties.If are the Laplace transform of the reactions and relative displacements, respectively, then where The free-motion equations in the Laplace domain can be obtained assembling the mass and the stiffness matrices associated with the structural configuration shown in Figure 7, resulting in where M = I 5 and The viscoelastic model is controlled by three parameters directly related to the damping, say ,   , and .In multipledof systems it seems logical to start a perturbation-based numerical method from the undamped problem.Therefore, checking the damping function, we note that the particular values   = 0 and  = 1 transform the nonlinear eigenvalue problem of (86) into the undamped linear one; that is, for these particular values, it is verified that G() ≡ 0. Among these parameters we chose  as our perturbation parameter.Initially, without previous information, there is no reason to take one or another parameter as long as it admits certain initial value associated with the undamped problem.However, the behavior of each parameter is different and depends on its roll within the mathematical viscoelastic model.A deeper research on the suitability of the different damping parameters of each damping model lies far from the objective of this paper, although it would be desirable as further research.
Considering then the storage parameter as variable, the damping function can be written as G(, ) (note that now our parameter is  ≡ ).Solving the undamped eigenvalue problem we can extract the natural frequencies and the normal modes {±  ,   } 5 =1 .Thus, the th modal damping function is Γ  (, ) =    G (, )   =  ()    Q  ≡  ()   .(88) After some straight operations the function Ψ  (, ) can be expressed as Moreover, considering the eigenvalues as -dependent functions, the first-and second-order approximations based on the Taylor series expansion are As aforementioned, the current method improves these approximations proposing the so-named ILS and IQS solutions whose expressions for this particular case are where the notation  has been used to denote the storage coefficient inside the integral.
In this example the results of the proposed method will be compared not only with the exact solution but also with another numerical method proposed by Adhikari and Pascual [18] with similar characteristics to the current: first, the eigenvalues are calculated by a noniterative approach; and, second, the method is valid for proportional or lightly nonproportional damping.Adhikari and Pascual's method assumes that the th eigenvalue can be written as   =  0 + , with  0 being certain point adequately chosen and  being an unknown value to be obtained.It is assumed that |/ 0 | ≪ 1; hence the solution is obtained expanding the characteristic modal equation D  () from (76) around  0 up to the first or up to the second order.Thus, the first-order approximation yields (using the same notation as [18]) while the second-order approximation yields where The lowest of the two roots (in absolute value) is chosen since a priori the correct value is not known.Adhikari and Pascual proposed as initial guess  0 = −    +   √1 −  2  , where   = lim  → 0 Γ  ()/2  , provided that this limit exists.However, when the viscoelastic function is based on fractional derivatives, this limit effectively is not a finite number.In fact, Adhikari and Pascual developed the method for Biot's model with multiple exponential kernels.This mathematical singularity can be avoided just taking  0 =   , that is, the undamped natural frequency, something that additionally helps us in the comparison with ILS and IQS approaches, since both of them also have the undamped eigenfrequency as initial point.
As shown in the mathematical results, the accuracy of the proposed approach depends directly on the distance between the damping parameter and its initial value.If this latter leads to the undamped problem, it is expected that the higher the parameter the more damped the system.For this reason, in the current example, two levels of damping will be considered separately: moderate damping (MD) and high damping (HD).Thus, the values of the damping parameters are chosen to represent both cases on the basis of the concept of loss factor peak,   , which can be used as a measure of the level of damping [42][43][44][45].According to the experimental evidence [44], values within the range 0.1 <   < 0.50 can be considered as moderate damping, while a loss factor 0.5 <   < 1 is representative of high damping materials, used for vibration control.According to this reasoning, we chose for the two cases the values of the parameters shown in Table 2.
The complex eigenvalues are computed using both methods (current and that of Adhikari-Pascual [18], in advance "AP") for each mode and for each level of damping.In order to arrange adequately the obtained results, the relative error (in percentage, %) with the exact solution of the modal characteristic equation is graphically represented in Figure 8. Moderate and high damping cases are shown in Figures 8(a) and 8(b), respectively.Within each plot, the type of solution (linear or quadratic) has been also included.Thus, the ILS and IQS have been calculated from (91) and (92) and the 1st-and 2nd-order solutions of AP method from (94) and (95), respectively.Both methods show very good agreement with the exact values for moderate damping.In fact, the relative error for every mode is always less than 0.01%.It seems logical to compare the ILS with the AP (1st-order approach) and IQS with the AP (2nd-order approach).In this regard, note that the proposed ILS method presents a relative error of the same order of magnitude as that of solutions based on the 2nd-order approximation and therefore can be considered as very good approximation.Additionally, the IQS exhibits the same level of accuracy as AP (2nd order), something that also occurs for high damping.Furthermore, for this particular case, ILS shows better results than AP (2nd order), but after several numerical examples with different damping parameters (not shown here) this does not always hold, although the order of magnitude of both approaches is the same.Another highlight is the loss of accuracy suffered by the AD (1st order), showing errors between 5% and 30% for high damping.Note that the relative error for all approaches increases with the mode; this is also in agreement with the fact that higher modes of proportional or lightly nonproportional systems become more damped than lower modes.Further research is currently being developed to cover with this method highly nonproportional systems.

Conclusions
A novel numerical method to compute the eigenvalues of viscoelastic oscillators is developed.The damping is introduced via a convolution integral with kernel functions which ensure that the dissipative forces depend on the history of the degreeof-freedom velocities.Each damping model is controlled by a kernel involving a set of parameters; the eigenvalues can be considered single-variable functions of some of these damping parameters.The application of the method requires finding only a particular value of the damping parameter and its corresponding eigenvalue.Under certain conditions, the eigenvalue can be expanded in its Taylor series and explicit expressions of the first and second derivatives are calculated.
The main contribution of this paper is the development of a numerical approach that improves the solution given by the Taylor series expansion.For that, it is shown that the characteristic equation can always be transformed into an ordinary differential equation with separated variables, so that the eigenvalue can be estimated by direct numerical integration.In addition, the order error of the proposed method is analyzed; while the th order Taylor series expansion presents O(ℎ +1 ), it is demonstrated that the current method presents O(ℎ +2 ).
To illustrate the results, the method is applied to the most commonly used damping models: viscous damping, nonviscous damping with exponential kernels, and fractional derivatives based models.For the first one, the method obtains approximated eigenvalues as analytical solutions of the proposed differential equation.For the second one, several forms of the damping parameters are presented, finding again analytical solutions for the Taylor expansions for all forms.The third one has been included as part of the viscoelastic constraints in a multiple degrees-of-freedom example.The numerical examples show that the lower the damping level is the better the accuracy of the proposed solutions is.The best results are always obtained when the parameter value is close to the considered initial point.Current research is being developed in two directions: (a) the implementation of the proposed method within systems with several damping functions of different nature and (b) the extension for nonproportionally damped multiple degreesof-freedom systems.

A. Proof of Property 1
The function Ψ(, ) is well defined, provided that the denominator         Γ (, ) +  (

B. Proof of Property 2
The function () is solution of the characteristic equation; therefore the following expression is verified:

Figure 1 :
Figure 1: Schematic representation of the viscoelastic oscillator.

Figure 2 :
Figure 2: Eigenvalues of viscous model versus .Continuous line for exact values and dotted one for Improved Solutions: square tick for Linear (ILS) and bullet for Quadratic (IQS).

Table 2 :
Viscoelastic model based on fractional derivatives.Values of the damping parameters for the two considered damping cases.