# Analysis of Nonviscous Oscillators Based on the Damping Model Perturbation

Authors: Mario Lázaro; César F. Casanova; Ignacio Ferrer; Pedro Martín
Journal: Shock and Vibration
 (2016)

---

## Abstract

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.


---

## Body

## 1. 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, sayu(t), can be modeled with a single mass, say m, attached to a fixed point by a viscoelastic constraint. Figure 1 shows the schematic configurationmass-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 R(t), is related with the displacement by(1)Rt=∫-∞tGt-τu˙τdτ+kut,where G(t) is the dissipative kernel or damping function, characteristic of the considered viscoelastic model, and k 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.Figure 1
Schematic representation of the viscoelastic oscillator.The motion equation can be deduced from the dynamic equilibrium:(2)mu¨+∫-∞tGt-τu˙τdτ+kut=Ftu0=u0,u˙0=v0,where u0 and v0 are the initial position and velocity of the mass and F(t) is the external, time-dependent, applied force.With reference to the mathematical form ofG(t), 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(t) 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 s. The form of this problem will depend on the nature of the viscoelastic function transform, named G(s)=L{G(t)}. 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 G(s) 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–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 G(s). 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 certaindesign 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 large-order 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, calledvalid 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.


## 2. Dynamics of Single Degree-of-Freedom Viscoelastic Systems

The characteristic equation of the eigenvalue problem can be obtained trying solutions with the formu(t)=u0est in the free-motion equation obtained from (2). Consider(3)ms2+sGs+k=0.Here, as mentioned, G(s) is the Laplace transform of the kernel function G(t). Note that the eigenvalue problem of proportional or lightly nonproportional MDOF systems of N degrees of freedom can be reduced to solve just N 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(4)Ds=s2+sΓs+ωn2=0.The damping function Γ(s)=G(s)/m and the natural frequency of the undamped system ωn=k/m have been introduced. Golla and Hughes [4] gave the necessary conditions that G(s) has to fulfil in order to define a real dissipative motion. Assuming that G(s) verifies these conditions, it can be assured [32, 33] that the characteristic equation (4) has 2+q eigenvalues, with q≥0. These can be arranged in a set with the form(5)s0,s0∗,s1,…,sq,where s0,s0∗ are a pair of complex conjugate roots which induce oscillatory motion. The rest {sj}j=1q are namednonviscous 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(6)ut=γ0a0t+γ0∗a0∗t+∑j=1qγjajt,where the function and coefficient are(7)ajt=∫0tesjt-τFτ-Gτu0dτ+mesjtv0+sju0,γj=∂Dsj∂s-1,0≤j≤q.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.


## 3. Perturbation Method of the Damping Model

### 3.1. Fundamentals of the Method

As explained in Section2, the eigenvalues are the roots of the characteristic equation given by (4). The key idea of the method is to consider that the function Γ=Γ(s,p1,…,pr) depends not only on the variable s but also on a set of damping parameters {p1,…,pr} which controls the viscoelastic behavior of the system.In the general case, any eigenvalueλ=λ(p1,…,pr) 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 p1. Denoting as θ=p1, 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 p2,p3,…,pr, obtaining the related initial eigenvalue λ0: (8)λ02+λ0Γλ0,θ0,p2,…,pr+ωn2=0⟶λ0=λ0θ0,p2,…,pr.Under these conditions, the parameter θ=p1 is namedvalid 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 (9)Dλ,θ=λ2+λΓλ,θ+ωn2=0defines an implicit function λ=λ(θ) around the point θ=θ0 if (10)∂Dλ0,θ0∂λ=Γλ0,θ0+λ02+∂Γλ0,θ0∂λ≠0.Under this hypothesis, it is also guaranteed that there exist (i) an open intervalI=θ0-h/2,θ0+h/2 with h>0 and (ii) an open neighborhood around the point λ0 in the complex domain B⊂C, such that the function λ(θ) will be unique and continuously differentiable in I and its graph verifies (11)θ,λθ:θ∈I=s,θ∈B×I:Ds,θ=0.Notice that the damping function has been represented directly as Γ(s,θ) in the implicit function theorem due to the fact that the parameters p2,…,pr 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 kth-order Taylor series of λ(θ) will be developed. Moreover, this approximation will be used in a further step to improve the solution.


### 3.2. 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 kth order can be obtained around the initial point θ0, provided that the derivatives evaluated at this point are available. Consider(12)λkθ=λ0+∑r=1kλrθ0r!θ-θ0rwith(13)λrθ0=drλdθrθ=θ0.The computation of λ(r)(θ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 (14)dDdθ=2λλ′+λ′Γ+λΓλ′λ′+Γθ′=0whence solving λ′ results in(15)λ′=-λΓθ′Γ+λ2+Γλ′,(16)Γλ′=∂Γ∂λ,Γθ′=∂Γ∂θ.Using the same procedure to compute the second derivative, the termλ′′ is extracted. After some simplifications, this second derivative results in (17)λ′′=-Rθ,λ,λ′Γ+λ2+Γλ′2,where (18)Rθ,λ,λ′=λ′22+2Γλ′+λΓλλ′′+2λ′Γθ′+λΓλθ′′+λΓθθ′′and now (19)Γλλ′′=∂2Γ∂λ2,Γλθ′′=∂2Γ∂λ∂θ,Γθθ′′=∂2Γ∂θ2.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.


### 3.3. 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(20)λ′=-λΨλ,θ,where the new function Ψ(λ,θ) is defined by (21)Ψλ,θ=Γθ′λ,θΓλ,θ+λ2+Γλ′λ,θ.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 kth-order solution of the Taylor series obtained in (12), namely, λk(θ), will be used to approximate the value of λ inside the function Ψ(λ,θ). So (22)Ψλ,θ≈Ψλkθ,θ≡ψkθ.Hence, the new function ψk(θ) 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 Γ(s,θ) with respect to s; a viscoelastic system with low viscoelasticity does not present large variations in s-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 Ψ(λ,θ)≈Ψ(λk(θ),θ) 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(23)λ′λ=-ψkθ,λθ0=λ0whose solution will be denoted by λ^k(θ). The subindex indicates the k-order Taylor series approximation. An explicit expression can be represented by (24)λ^kθ=λ0e-Zθ,where the exponential function is defined by(25)Zθ=∫θ0θψkϑdϑ.When the functionλ^k(θ) is obtained from the first-order approximation of the Taylor series expansion k=1, the solution will be namedImproved Linear Solution (ILS). If the second-order approximation k=2 is used, the name will beImproved Quadratic Solution (IQS). The functions Ψ(s,θ) and ψk(θ) will be called exact and approximatedintegrand functions, respectively. The next subsection is aimed at the analysis of the truncation error.


## 3.1. Fundamentals of the Method

As explained in Section2, the eigenvalues are the roots of the characteristic equation given by (4). The key idea of the method is to consider that the function Γ=Γ(s,p1,…,pr) depends not only on the variable s but also on a set of damping parameters {p1,…,pr} which controls the viscoelastic behavior of the system.In the general case, any eigenvalueλ=λ(p1,…,pr) 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 p1. Denoting as θ=p1, 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 p2,p3,…,pr, obtaining the related initial eigenvalue λ0: (8)λ02+λ0Γλ0,θ0,p2,…,pr+ωn2=0⟶λ0=λ0θ0,p2,…,pr.Under these conditions, the parameter θ=p1 is namedvalid 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 (9)Dλ,θ=λ2+λΓλ,θ+ωn2=0defines an implicit function λ=λ(θ) around the point θ=θ0 if (10)∂Dλ0,θ0∂λ=Γλ0,θ0+λ02+∂Γλ0,θ0∂λ≠0.Under this hypothesis, it is also guaranteed that there exist (i) an open intervalI=θ0-h/2,θ0+h/2 with h>0 and (ii) an open neighborhood around the point λ0 in the complex domain B⊂C, such that the function λ(θ) will be unique and continuously differentiable in I and its graph verifies (11)θ,λθ:θ∈I=s,θ∈B×I:Ds,θ=0.Notice that the damping function has been represented directly as Γ(s,θ) in the implicit function theorem due to the fact that the parameters p2,…,pr 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 kth-order Taylor series of λ(θ) will be developed. Moreover, this approximation will be used in a further step to improve the solution.


## 3.2. 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 kth order can be obtained around the initial point θ0, provided that the derivatives evaluated at this point are available. Consider(12)λkθ=λ0+∑r=1kλrθ0r!θ-θ0rwith(13)λrθ0=drλdθrθ=θ0.The computation of λ(r)(θ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 (14)dDdθ=2λλ′+λ′Γ+λΓλ′λ′+Γθ′=0whence solving λ′ results in(15)λ′=-λΓθ′Γ+λ2+Γλ′,(16)Γλ′=∂Γ∂λ,Γθ′=∂Γ∂θ.Using the same procedure to compute the second derivative, the termλ′′ is extracted. After some simplifications, this second derivative results in (17)λ′′=-Rθ,λ,λ′Γ+λ2+Γλ′2,where (18)Rθ,λ,λ′=λ′22+2Γλ′+λΓλλ′′+2λ′Γθ′+λΓλθ′′+λΓθθ′′and now (19)Γλλ′′=∂2Γ∂λ2,Γλθ′′=∂2Γ∂λ∂θ,Γθθ′′=∂2Γ∂θ2.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.


## 3.3. 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(20)λ′=-λΨλ,θ,where the new function Ψ(λ,θ) is defined by (21)Ψλ,θ=Γθ′λ,θΓλ,θ+λ2+Γλ′λ,θ.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 kth-order solution of the Taylor series obtained in (12), namely, λk(θ), will be used to approximate the value of λ inside the function Ψ(λ,θ). So (22)Ψλ,θ≈Ψλkθ,θ≡ψkθ.Hence, the new function ψk(θ) 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 Γ(s,θ) with respect to s; a viscoelastic system with low viscoelasticity does not present large variations in s-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 Ψ(λ,θ)≈Ψ(λk(θ),θ) 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(23)λ′λ=-ψkθ,λθ0=λ0whose solution will be denoted by λ^k(θ). The subindex indicates the k-order Taylor series approximation. An explicit expression can be represented by (24)λ^kθ=λ0e-Zθ,where the exponential function is defined by(25)Zθ=∫θ0θψkϑdϑ.When the functionλ^k(θ) is obtained from the first-order approximation of the Taylor series expansion k=1, the solution will be namedImproved Linear Solution (ILS). If the second-order approximation k=2 is used, the name will beImproved Quadratic Solution (IQS). The functions Ψ(s,θ) and ψk(θ) will be called exact and approximatedintegrand functions, respectively. The next subsection is aimed at the analysis of the truncation error.


## 4. Analysis of the Error

Section3 was focused on the methodology exposition to compute the function λ^k(θ) 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 ψk(θ) is not important. In this subsection more quantitative information is provided, finding a bound of the error |λ^k(θ)-λ(θ)|. The necessary conditions to calculate this bound are directly related with the properties of the viscoelastic function G(s,θ)=mΓ(s,θ) and will be given by three hypotheses:(H1)
If(λ0,θ0) is the initial point verifying (8), then Γ(λ0,θ0)+λ0[2+Γλ′(λ0,θ0)]≠0.(H2)
LetI and B be the neighborhoods of θ0 and λ0, respectively; the existence of these sets is assured by the implicit function theorem. Then G(s,θ) is a continuous and differentiable function with respect to the variables s,θ up to order k+1 in the set B×I. In addition, G(s,θ) and its partial derivatives with respect to s,θ up to order k+1 are bounded in the set B.(H3)
The functionΨ(s,θ) defined by (21) is Lipschitz continuous with respect to the first variable s in the set B, with constant q>0: (26)Ψs1,θ-Ψs2,θ≤qs1-s2,∀s1,s2∈B.The hypothesis (H1) allows assuring the existence of the functionλ(θ):I→B in the interval I=θ0-h/2,θ0+h/2. One of the theses of the implicit function theorem is the existence of the set B⊂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.Property 1.
IfΨ(s,θ) is a function from (21), then a real and positive number exists, p>0, such that |Ψ(s,θ)|≤p∀(s,θ)∈B×I.Property 2.
Letλ(θ) be the eigenvalue as function of the parameter θ∈I. Then for any integer k≥0 there exists a real and positive number, ρk>0, such that |λ(k)(θ)|≤ρk∀θ∈I.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 λ^k(θ) be the approximation calculated from (22) to (25) with k≥1. Assuming that hypotheses (H1), (H2), and (H3) are satisfied, there exist two positive real numbers, ρk>0 and p>0, such that (27)λ^kθ-λθ≤qλ0ρkk+1!hk+2eph,∀θ∈I.Proof.
From (20), the exact eigenvalue λ(θ) verifies the relation (28)λθ=λ0e-Yθ,where (29)Yθ=∫θ0θΨλϑ,ϑdϑ.
For demonstration purposes, it is necessary to bound the difference norm|Zn(θ)-Yn(θ)| for any natural number n≥1.
First, using (H3) and the subtractionλ^k(θ)-λ(θ) in Taylor series, the previous norm for n=1 is(30)Zθ-Yθ≤∫θ0θΨλkϑ,ϑ-Ψλϑ,ϑdϑ≤∫θ0θqλkϑ-λϑdϑ≤q∫θ0θλk+1ξk+1!ξ-θ0k+1dϑ≤qλk+1ξk+1!hk+2≤qρkk+1!hk+2,where ξ∈I defines the truncation residual of the Taylor series up to the order k. For the last inequality Property 2 has been used. In general, a bound of the functions Z(θ) and Y(θ) can be calculated using Property 1. From the definitions, one directly obtains (31)Zθ≤∫θ0θΨλkϑ,ϑdϑ≤ph,Yθ≤∫θ0θΨλϑ,ϑdϑ≤ph.
Second, using (30) and (31), for n≥1,(32)Znθ-Ynθ=Zθ-Yθ∑j=0nZn-jθYj-1θ≤Zθ-Yθ∑j=1nZθn-jYθj-1≤Zθ-Yθ∑j=1nphn-jphj-1=Zθ-Yθnphn-1.
Now the error can be calculated; using the definitions of both functions,(33)λ^kθ-λθ=λ0∑n=1∞-1nn!Znθ-Ynθ≤λ0∑n=1∞1n!Znθ-Ynθ≤λ0Zθ-Yθ∑n=1∞nphn-1n!≤λ0Zθ-Yθ∑n=1∞phn-1n-1!≤λ0Zθ-Yθeph≤qλ0ρkk+1!hk+2eph.The theorem states that the approximation order of the improved solution isO(hk+2), while approximating with Taylor series up to the kth term gives O(hk+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.


## 5. Numerical Examples

### 5.1. The Viscous Model

The viscous model is a special case of viscoelastic behavior whose kernel is proportional to the Dirac delta,G(t)=cvδ(t), where cv is the viscous damping coefficient. The Laplace transform of the kernel function becomes the s-independent function G(s)≡cv. In general, instead of cv the damping ratio ζ is commonly used. Both are related by the well-known expression cv=2mωnζ, and consequently the mathematical model depends only on the single parameter θ≡ζ. The mass-normalized characteristic equation (3) can be rewritten as(34)s2+2sζωn+ωn2=0,where now Γ(s,ζ)=2ζωn. The classical exact solution of (34) is the pair of complex numbers (35)λζ=-ωnζ±iωn1-ζ2.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=±iωn. Using now the eigenvalue with positive imaginary part, the derivatives of λ(ζ) can be deduced from (15) and (17): (36)λ′ζ=-ωnλζωnζ+λζ⟶λ′0=-ωn,λ′′ζ=-λ′ζ2ωn+λ′ζωnζ+λζ⟶λ′′0=-iωn.The Taylor series expansion of the functionλ(ζ) up to the second order can be obtained as (37)λ2ζ=iωn-ωnζ-iωnζ22≡λ1ζ-iωnζ22.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 (k=1) and IQS (k=2) can be computed solving the following differential equation: (38)λ′=-λωnωnζ+λkζ,λ0=iωn.Substituting in (38) the first order approximation of the root λ1(ζ) from (37), the ordinary differential equation leads to (39)λ′=iλ,λ0=iωnand a closed expression for ILS is obtained: (40)λ^1ζ=iωneiζ.It is also interesting to deduce the improved solution corresponding to the quadratic approximation λ2(ζ) from (37); the differential equation results in (41)λ′=λi1-ζ2/2,λ0=iωnand the explicit expression of IQS can also be easily integrated, resulting in (42)λ^2ζ=iωn2+ζ2-ζi/2.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 ωn=1 rad/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.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).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): (43)λ^1ζ=λ^2ζ=λζ=ωn,∀ζ≥0which coincides with the absolute value of the exact root for any damping ratio.


### 5.2. 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 (44)Gt=cvμe-μt,μ,t≥0,where μ is the relaxation parameter, also callednonviscous parameter. Notice that for high values of μ the viscoelastic model tends to the viscous model, mathematically expressed by limμ→∞⁡G(t)=cvδ(t). Therefore, the nonviscous behavior will be controlled by μ, so that when μ≫ωn the damping model is viscous with coefficient cv [34, 35]. Otherwise, when μ≪ωn the behavior is highly nonviscous. The characteristic equation in the Laplace domain can be written as (45)ms2+sGs+k=0,where G(s)=cvμ/(s+μ). The coefficient of the limit viscous model can again be expressed in terms of the damping ratio cv=2mωnζ. Under these conditions, the viscoelastic model can be written as function of the two parameters, ζ and μ, as follows: (46)λ2+2ωnζμλ+μλ+ωn2=0.Introducing a new parameterτ=1/μ with μ>0, the previous characteristic equation can be rewritten as(47)λ2+2ωnζ1τλ+1λ+ωn2=0.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 N exponential kernels in Section 5.3.Table 1
Biot’s exponential model. Damping parameters, considered eigenvalue functions, and initial values of the proposed method.
Parameter,θ
Functionλθ
Initial value,θ0
Initial eigenvalue,λ0=λθ0
Section
D
(
λ
,
θ
)
Results
θ
=
μ
λ
=
λ
μ
μ
0
=
0
λ
0
=
±
i
ω
n
Section5.2 (Case 1)
Equation (48)
Figure3
θ
=
τ
λ
=
λ
τ
τ
0
=
0
λ
0
=
-
ω
n
ζ
±
i
ω
n
1
-
ζ
2
Section5.2 (Case 2)
Equation (53)
Figure4
θ
=
ζ
λ
=
λ
ζ
ζ
0
=
0
λ
0
=
±
i
ω
n
Section5.3
Equation (60)
Figure5Case 1 (θ=μ).
The functionλ(μ) is implicitly defined by(48)Dλ,μ=λ2+2ωnζμλ+μλ+ωn2=0.As an initial point, any root of the equation D(λ,0)=λ2+ωn2=0 can be taken. Without loss of generality we will use the positive imaginary part λ0=iω. From (15) and (17) the values of λ′(0),λ′′(0) can be computed, obtaining after some operations the expressions(49)λ′0=iζ,λ′′0=-ζ2+ζωn.Hence, the second-order Taylor polynomial is (50)λ2μ=iωn+iζμ-ζ2+ζ2ωnμ2,μ≥0.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: (51)ψkμ=λkμζωnλkμλkμ+μ2+ζμ2ωn,where λk(μ) is defined by (12). As described before, the quality of the approximated eigenvalue function λ^k(μ) depends on the accuracy of the assumption over the integrand function Ψλ(μ),μ≈ψkμ 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 ωn=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 Ψ(λ,μ),(52)ϵkμ=Ψλμ,μ-ψkμΨλμ,μ,is shown. Notice that the exact eigenvalue λ(μ) from (46) is required for the error computation.Figure 3
Caseθ=μ. Relative error ϵk(μ)=Ψ-ψk/Ψ of the integrand function (a). Eigenvalues versus parameter μ (b and c). Continuous line for exact values and dotted one for Improved Solutions: square tick for Linear (ILS) and bullet for Quadratic (IQS).
(a)
(b)
(c)Figure3 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(53)Dλ,τ=λ2+2ωnζ1τλ+1λ+ωn2=0.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:(54)λ′0=ζλ02ωnλ0+ζωn,λ′′0=-ζλ03ωn2λ02+ζλ0ωn-2ζ2ωn2λ0+ζωn3.With these expressions, the second-order Taylor series expansion of λ(τ) is (55)λ2τ=λ0+ζλ02ωnλ0+ζωnτ-2λ02+ζλ0ωn-2ζ2ωn22λ0+ζωn3ζλ03ωnτ2defined in the range ζ,τ≥0. Now, ILS and IQS can be calculated with the methodology described in Section 3.3. First, in this particular case, (56)ψkτ=-λkτζωnλkτ1+τλkτ2+ζωn.Second, for each value of the parameter τ, the eigenvalue λ^k(τ)=λ0e-Z(τ) can be estimated by numerical quadrature of the integral:(57)Zτ=∫0τψktdt.
Figure4 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.Figure 4
Caseθ=τ=1/μ. Relative error ϵk(τ)=Ψ-ψk/Ψ of the integrand function (a). Eigenvalues versus parameter μ (b and c). Continuous line for exact values and dotted one for Improved Solutions: square tick for Linear (ILS) and bullet for Quadratic (IQS).
(a)
(b)
(c)


### 5.3. 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(58)Gt=cvN∑j=1Nμje-μjt,μj,t≥0,1≤j≤N,where μj are the relaxation parameters and cv=2mωnζ is, as before, the damping coefficient of the viscous model when μj→∞. The Laplace transform of (58) can be expressed as (59)Gs,ζ,μ1,…,μN=2mωnζN∑j=1Nμjs+μj.Directly, the mass-normalized characteristic equation takes the form (60)λ2+2λωnζN∑j=1Nμjλ+μj+ωn2=0.Note that the previous equation can be transformed into a polynomial of order 2+N. The set formed by 2+N eigenvalues can be separated, on one side, in a subset with the two complex conjugate eigenvalues and, on the other, in another with N nonviscous eigenvalues that in general are real and negative. The parameters μj cannot be used with the proposed method because analytical solutions at the initial point for N 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=±iωn. 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, (61)λ′0=-ωnN∑j=1Nμjλ0+μj,λ′′0=λ′0λ′0λ0+2N∑j=1Nωnμjλ0+μj2.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(62)ψkζ=ζ+λkζ1-ζωn/N∑j=1Nμj/λkζ+μj2ωn/N∑j=1Nμj/λkζ+μj-1.The method is implemented for a numerical example withN=5 kernels. The relaxation parameters are taken as μj={2,4,6,10,15} rad/s. Three different natural frequencies ωn={0.2,2,20} rad/s have been considered to validate the results when the system is near-viscous μj≫ωn or, otherwise, nonviscous μj≈ωn.Figure5 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 ωn=0.20 rad/s (viscous behavior), where the slope of the imaginary part becomes almost vertical.Figure 5
Caseθ=ζ for Biot’s multiexponential model μj={2,4,6,10,15} rad/s. Relative error ϵk(ζ)=Ψ-ψk/Ψ of the integrand function (a). Eigenvalues versus parameter ζ (b and c). Continuous line for exact values and dotted one for Improved Solutions: square tick for Linear (ILS) and bullet for Quadratic (IQS).
(a)
(b)
(c)Real Eigenvalues. It is expected that the real eigenvalues of (59) will be relatively close to the parameters μj in lightly damped systems, as described by Adhikari and Pascual [19]. Naming λ=λ(ζ) to the nonviscous eigenvalue associated with a particular parameter μr, then λ0=limζ→0⁡λ(ζ)=-μr. Indeed, multiplying (60) by λ+μr results in(63)λζ+μrλ2ζ+ωn2+2λζωnζNμr+2λζωnζN∑j=1j≠rNλζ+μrλζ+μjμj=0.Taking limits ζ→0 in the previous expression,(64)λ0+μrλ02+ωn2=0.Since λ0∈R for the nonviscous eigenvalues, λ02+ωn2>0 and consequently λ0+μr=0. Therefore, the nonviscous eigenvalue associated with μj can be expanded around the point λ0=λ(0)=-μj. 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 s=-μj (see (59)), the derivatives from (15) and (17) cannot be calculated. To avoid this problem, (63) with r=1 will be used. Thus, λ′(0) and λ′′(0) can be calculated taking derivatives with respect to ζ in the expression:(65)λζ+μ1λ2ζ+ωn2+2λζωnζNμ1+∑j=2Nλζ+μ1λζ+μjμj=0. Evaluating (65) at ζ=0, and after some operations, the following results are obtained: (66)λ′0=2ωnμ12ωn2+μ12,λ′′0=8ωn2μ13ωn2+μ122μ12-ωn2ωn2+μ12+∑j=2Nμjμj-μ1.The computation of the improved solutions from the differential equation (20) requires that G(s) is well defined at the initial point. Since λ0=-μj is a pole of G(s) and of ∂G(s)/∂s, the Taylor series up to the first order and second order will be used to estimate the eigenvalues.Figure6 shows results for a numerical example with ωn=3 rad/s and Biot’s model with N=3 kernels and μj={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,λIV0,…, from (65) can be extracted and new terms can be added.Figure 6
Exact solution and Taylor series approximation (quadratic and linear) of the nonviscous eigenvalues versus damping ratioζ. Biot’s model with N=3 kernels and nonviscous parameters μj={1,5,10} rad/s.


### 5.4. Multiple-Degree-of-Freedom Systems

In this last point, the application of the proposed method for multiple degrees-of-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 considern degrees-of-freedom (dof) vibrating system. The dofs’ time-domain response is represented by a vector u(t)∈Rn. With help of the Finite Element Method, the mass and the stiffness matrices of the system, denoted, respectively, by M,K∈Rn×n, 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 fdt∈Rn are proportional to the history of the dofs’ velocities via kernel hereditary functions. These functions are the entries of a matrix denoted by G(t)∈Rn×n, also assumed symmetric. Consider(67)fdt=∫-∞tGt-τu˙τdτ.Assuming small displacements, the motion equations linearly relate the dynamic balance among inertiaMu¨, elastic Ku, damping fd(t), and external forces fe(t). The general form of the system of integrodiferential equations can be written by (68)Mu¨+∫-∞tGt-τu˙τdτ+Ku=fet,u˙0=v0,u0=u0.Under these conditions, the modes of the system can be obtained as the nontrivial solutions of the free-motion problem obtained considering fe(t)=v0=u0=0 in (68). Thus, checking functions of the form u(t)=uest, we obtain (69)s2M+sGs+Ku≡Dsu=0,where G(s)=L{G(t)}∈Cn×n represents the Laplace transform of the viscoelastic kernels and D(s):C→Cn×n is thedynamic 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 (70)det⁡Ds=0.The complete set of eigenmodes of this problem present 2n+q distinct values distributed as(71)λ1,…,λn,λ1∗,…,λn∗,σ1,…,σq,where λj,λj∗∈C, 1 ≤j≤n, are n conjugate complex pairs, corresponding to the modes with oscillatory nature. The rest, σj∈R-, 1≤j≤q, are negative real numbers that represent overcritically damped modes callednonviscous 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 uj,uj∗∈Cn,1≤j≤n to the complex eigenvectors associated with λj,λj∗, and aj∈Rn,1≤j≤q to the eigenvector associated with the real eigenvalue σj.In the analysis and search of solutions of the eigenproblem, the proportionality ofn-dof system, that is, the modal decoupling capability, becomes of special importance. As known, the undamped problem obtained from G(t)≡0 can always be diagonalized in the modal space of matrices M and K. Denoting by ϕj∈Rn to the jth mass-homogenized real mode (in column) and by ωj to its associated natural frequency, the classical orthogonal relations can be written as (72)ϕjTMϕk=δjk,ϕjTKϕk=ωj2δjk,where δjk is the Kronecker delta. The modal matrix Φ groups in columns n undamped modes and is used to change the dof’s to the modal coordinates by u=Φq. Hence, (69) can be expressed as (73)s2IN+sΓs+Λq=0,where Λ=ΦTKΦ=diag⁡ωj2 is the diagonal squared natural frequency matrix and Γ(s)=ΦTG(s)Φ 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 calledproportional, since the necessary and sufficient conditions for the modal decoupling are directly related with proportional relationships between the damping matrix G(s) and the dynamic matrices M,K [31]. In some cases of nonproportionality, the matrix Γ(λj) is diagonally dominant but not purely diagonal, fact which is equivalent to assume as true: (74)∑l=1l≠knΓklλj<Γkkλj,∀1≤k,j≤n.The system is then said to belightly 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 (75)det⁡s2In+sΓs+Λ≈∏j=1ns2+sΓjjs+ωj2.Hence, the set of eigenvalues can be obtained from the following n decoupled equations: (76)Djs=s2+sΓjjs+ωj2=0,1≤j≤n,where Γjj(s)=ϕjTG(s)ϕj. Equation (76) represents the jth 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 Γjj(s,θ)=ϕjTG(s,θ)ϕj to refer to the jth modal damping function. For certain particular value of this parameter, say θ=θ0, (76) admits a closed-form analytic solution, denoted by λj(θ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 (77)λj1θ=λjθ0+λj′θ0θ-θ0,(78)λj2θ=λjθ0+λj′θ0θ-θ0+λj′′θ0θ-θ022,where the subscripts (·)j1 and (·)j2 denote the first- and second-order approximations of the jth eigenvalue, respectively. Expressions for λj′ and λj′′ 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 j: (79)Ψjs,θ=∂Γjj/∂θΓjjs,θ+s2+∂Γjj/∂sas a generalization of that of (21). It has been proved that theImproved Linear and Quadratic Solutions (ILS and IQS, see (24)) defined as (80)ILS:λ^j1θ=λjθ0exp⁡-∫θ0θΨjλj1ϑ,ϑdϑ,(81)IQS:λ^j2θ=λjθ0exp⁡-∫θ0θΨjλj2ϑ,ϑdϑ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 Figure7. 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 j=0,6 of fixed boundaries are related with zero dof’s so that u0=u˙0=u6=u˙6=0. Thus, the constitutive equations can be written as(82)Rj-1,j+TrαdαRj-1,jdtα=kΔuj+cTrαdαΔujdtα,j=1,3,4,5,where Δuj=uj-uj-1 and c,α, and Tr 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 c>1,0<α<1,Tr>0. The coefficient k, is the linear, static rigidity.Figure 7
Multiple degrees-of-freedom example: lumped-mass dynamical system with viscoelastic links based on the fractional derivatives.The kernel functionG(t) 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 G(s) can easily be calculated by simply applying to the fractional derivatives of (82) the Laplace transform and using its properties. If (83)Rj-1,js=LRj-1,jt,Δu^js=u^js-u^j-1s=LΔu^jtare the Laplace transform of the reactions and relative displacements, respectively, then (84)Rjs=k1+cTrsα1+TrsαΔu^js≡k+sGsΔu^js,where(85)Gs=ksc-1Trsα1+Trsα.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 Figure7, resulting in(86)s2M+sGs+Ku^s=0,where M=mI5 and (87)K=k2-1000-12-1000-12-1000-12-1000-12,Gs=Gs1000001-1000-12-1000-12-1000-11≡GsQ.The viscoelastic model is controlled by three parameters directly related to the damping, sayc, Tr, and α. In multiple-dof 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 Tr=0 and c=1 transform the nonlinear eigenvalue problem of (86) into the undamped linear one; that is, for these particular values, it is verified that G(s)≡0. Among these parameters we chose c 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 asG(s,c) (note that now our parameter is θ≡c). Solving the undamped eigenvalue problem we can extract the natural frequencies and the normal modes {±ωj,ϕj}j=15. Thus, the jth modal damping function is(88)Γjjs,c=ϕjTGs,cϕj=GsϕjTQϕj≡GsQj.After some straight operations the function Ψj(s,c) can be expressed as (89)Ψjs,c=∂Γjj/∂θΓjjs,θ+s2+∂Γjj/∂s=1+sTrααc-1+2s21+sTrα2/kQjsTrα.Moreover, considering the eigenvalues as c-dependent functions, the first- and second-order approximations based on the Taylor series expansion are (90)λj1c=λj1+λj′1c-1=iωj-c-1kQjTriωjTrα-121+iωjTrα,λj2c=λj1+λj′1c-1+λj′′1c-122=λj1c-ic-12k2Qj2iωjTr2α1-2α+iωjTrα8ωj31+iωjTrα3.As aforementioned, the current method improves these approximations proposing the so-named ILS and IQS solutions whose expressions for this particular case are (91)ILS:λ^j1c=iωjexp⁡-∫γ=1cΨjλj1γ,γdγ,(92)IQS:λ^j2c=iωjexp⁡-∫γ=1cΨjλj2γ,γdγ,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 jth eigenvalue can be written as sj=s0j+δ, with s0j being certain point adequately chosen and δ being an unknown value to be obtained. It is assumed that δ/s0j≪1; hence the solution is obtained expanding the characteristic modal equation Dj(s) from (76) around s0j up to the first or up to the second order. Thus, the first-order approximation yields (using the same notation as [18]) (93)sj1AP≈s0j+δ1=s0j-Djs0j∂Djs0j/∂s,while the second-order approximation yields(94)sj2AP≈s0j+δ2=s0j+-B±B2-4AC2A,where (95)A=12∂2Djs0j∂s2,B=∂Djs0j∂s,C=Djs0j.The lowest of the two roots (in absolute value) is chosen sincea priori the correct value is not known. Adhikari and Pascual proposed as initial guess s0j=-ξjωj+iωj1-ξj2, where ξj=lims→0⁡Γjj(s)/2ωj, 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 s0j=iωj, 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 thedistance 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 ofloss factor peak, ηm, which can be used as a measure of the level of damping [42–45]. According to the experimental evidence [44], values within the range 0.1<ηm<0.50 can be considered as moderate damping, while a loss factor 0.5<ηm<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.Table 2
Viscoelastic model based on fractional derivatives. Values of the damping parameters for the two considered damping cases.
Damping level
c
α
T
r (sec.)
η
m
Moderate damping (MD)
10
0,50
1
0
-
6
0,45
High damping (HD)
500
0,50
1
0
-
6
0,90The 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.Figure 8
Relative error in the computation of eigenvalues: proposed (solid lines) and Adhikari-Pascual [18] (dashed lines). (a) Moderate damping. (b) High damping.
(a)
(b)


## 5.1. The Viscous Model

The viscous model is a special case of viscoelastic behavior whose kernel is proportional to the Dirac delta,G(t)=cvδ(t), where cv is the viscous damping coefficient. The Laplace transform of the kernel function becomes the s-independent function G(s)≡cv. In general, instead of cv the damping ratio ζ is commonly used. Both are related by the well-known expression cv=2mωnζ, and consequently the mathematical model depends only on the single parameter θ≡ζ. The mass-normalized characteristic equation (3) can be rewritten as(34)s2+2sζωn+ωn2=0,where now Γ(s,ζ)=2ζωn. The classical exact solution of (34) is the pair of complex numbers (35)λζ=-ωnζ±iωn1-ζ2.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=±iωn. Using now the eigenvalue with positive imaginary part, the derivatives of λ(ζ) can be deduced from (15) and (17): (36)λ′ζ=-ωnλζωnζ+λζ⟶λ′0=-ωn,λ′′ζ=-λ′ζ2ωn+λ′ζωnζ+λζ⟶λ′′0=-iωn.The Taylor series expansion of the functionλ(ζ) up to the second order can be obtained as (37)λ2ζ=iωn-ωnζ-iωnζ22≡λ1ζ-iωnζ22.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 (k=1) and IQS (k=2) can be computed solving the following differential equation: (38)λ′=-λωnωnζ+λkζ,λ0=iωn.Substituting in (38) the first order approximation of the root λ1(ζ) from (37), the ordinary differential equation leads to (39)λ′=iλ,λ0=iωnand a closed expression for ILS is obtained: (40)λ^1ζ=iωneiζ.It is also interesting to deduce the improved solution corresponding to the quadratic approximation λ2(ζ) from (37); the differential equation results in (41)λ′=λi1-ζ2/2,λ0=iωnand the explicit expression of IQS can also be easily integrated, resulting in (42)λ^2ζ=iωn2+ζ2-ζi/2.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 ωn=1 rad/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.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).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): (43)λ^1ζ=λ^2ζ=λζ=ωn,∀ζ≥0which coincides with the absolute value of the exact root for any damping ratio.


## 5.2. 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 (44)Gt=cvμe-μt,μ,t≥0,where μ is the relaxation parameter, also callednonviscous parameter. Notice that for high values of μ the viscoelastic model tends to the viscous model, mathematically expressed by limμ→∞⁡G(t)=cvδ(t). Therefore, the nonviscous behavior will be controlled by μ, so that when μ≫ωn the damping model is viscous with coefficient cv [34, 35]. Otherwise, when μ≪ωn the behavior is highly nonviscous. The characteristic equation in the Laplace domain can be written as (45)ms2+sGs+k=0,where G(s)=cvμ/(s+μ). The coefficient of the limit viscous model can again be expressed in terms of the damping ratio cv=2mωnζ. Under these conditions, the viscoelastic model can be written as function of the two parameters, ζ and μ, as follows: (46)λ2+2ωnζμλ+μλ+ωn2=0.Introducing a new parameterτ=1/μ with μ>0, the previous characteristic equation can be rewritten as(47)λ2+2ωnζ1τλ+1λ+ωn2=0.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 N exponential kernels in Section 5.3.Table 1
Biot’s exponential model. Damping parameters, considered eigenvalue functions, and initial values of the proposed method.
Parameter,θ
Functionλθ
Initial value,θ0
Initial eigenvalue,λ0=λθ0
Section
D
(
λ
,
θ
)
Results
θ
=
μ
λ
=
λ
μ
μ
0
=
0
λ
0
=
±
i
ω
n
Section5.2 (Case 1)
Equation (48)
Figure3
θ
=
τ
λ
=
λ
τ
τ
0
=
0
λ
0
=
-
ω
n
ζ
±
i
ω
n
1
-
ζ
2
Section5.2 (Case 2)
Equation (53)
Figure4
θ
=
ζ
λ
=
λ
ζ
ζ
0
=
0
λ
0
=
±
i
ω
n
Section5.3
Equation (60)
Figure5Case 1 (θ=μ).
The functionλ(μ) is implicitly defined by(48)Dλ,μ=λ2+2ωnζμλ+μλ+ωn2=0.As an initial point, any root of the equation D(λ,0)=λ2+ωn2=0 can be taken. Without loss of generality we will use the positive imaginary part λ0=iω. From (15) and (17) the values of λ′(0),λ′′(0) can be computed, obtaining after some operations the expressions(49)λ′0=iζ,λ′′0=-ζ2+ζωn.Hence, the second-order Taylor polynomial is (50)λ2μ=iωn+iζμ-ζ2+ζ2ωnμ2,μ≥0.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: (51)ψkμ=λkμζωnλkμλkμ+μ2+ζμ2ωn,where λk(μ) is defined by (12). As described before, the quality of the approximated eigenvalue function λ^k(μ) depends on the accuracy of the assumption over the integrand function Ψλ(μ),μ≈ψkμ 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 ωn=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 Ψ(λ,μ),(52)ϵkμ=Ψλμ,μ-ψkμΨλμ,μ,is shown. Notice that the exact eigenvalue λ(μ) from (46) is required for the error computation.Figure 3
Caseθ=μ. Relative error ϵk(μ)=Ψ-ψk/Ψ of the integrand function (a). Eigenvalues versus parameter μ (b and c). Continuous line for exact values and dotted one for Improved Solutions: square tick for Linear (ILS) and bullet for Quadratic (IQS).
(a)
(b)
(c)Figure3 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(53)Dλ,τ=λ2+2ωnζ1τλ+1λ+ωn2=0.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:(54)λ′0=ζλ02ωnλ0+ζωn,λ′′0=-ζλ03ωn2λ02+ζλ0ωn-2ζ2ωn2λ0+ζωn3.With these expressions, the second-order Taylor series expansion of λ(τ) is (55)λ2τ=λ0+ζλ02ωnλ0+ζωnτ-2λ02+ζλ0ωn-2ζ2ωn22λ0+ζωn3ζλ03ωnτ2defined in the range ζ,τ≥0. Now, ILS and IQS can be calculated with the methodology described in Section 3.3. First, in this particular case, (56)ψkτ=-λkτζωnλkτ1+τλkτ2+ζωn.Second, for each value of the parameter τ, the eigenvalue λ^k(τ)=λ0e-Z(τ) can be estimated by numerical quadrature of the integral:(57)Zτ=∫0τψktdt.
Figure4 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.Figure 4
Caseθ=τ=1/μ. Relative error ϵk(τ)=Ψ-ψk/Ψ of the integrand function (a). Eigenvalues versus parameter μ (b and c). Continuous line for exact values and dotted one for Improved Solutions: square tick for Linear (ILS) and bullet for Quadratic (IQS).
(a)
(b)
(c)


## 5.3. 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(58)Gt=cvN∑j=1Nμje-μjt,μj,t≥0,1≤j≤N,where μj are the relaxation parameters and cv=2mωnζ is, as before, the damping coefficient of the viscous model when μj→∞. The Laplace transform of (58) can be expressed as (59)Gs,ζ,μ1,…,μN=2mωnζN∑j=1Nμjs+μj.Directly, the mass-normalized characteristic equation takes the form (60)λ2+2λωnζN∑j=1Nμjλ+μj+ωn2=0.Note that the previous equation can be transformed into a polynomial of order 2+N. The set formed by 2+N eigenvalues can be separated, on one side, in a subset with the two complex conjugate eigenvalues and, on the other, in another with N nonviscous eigenvalues that in general are real and negative. The parameters μj cannot be used with the proposed method because analytical solutions at the initial point for N 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=±iωn. 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, (61)λ′0=-ωnN∑j=1Nμjλ0+μj,λ′′0=λ′0λ′0λ0+2N∑j=1Nωnμjλ0+μj2.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(62)ψkζ=ζ+λkζ1-ζωn/N∑j=1Nμj/λkζ+μj2ωn/N∑j=1Nμj/λkζ+μj-1.The method is implemented for a numerical example withN=5 kernels. The relaxation parameters are taken as μj={2,4,6,10,15} rad/s. Three different natural frequencies ωn={0.2,2,20} rad/s have been considered to validate the results when the system is near-viscous μj≫ωn or, otherwise, nonviscous μj≈ωn.Figure5 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 ωn=0.20 rad/s (viscous behavior), where the slope of the imaginary part becomes almost vertical.Figure 5
Caseθ=ζ for Biot’s multiexponential model μj={2,4,6,10,15} rad/s. Relative error ϵk(ζ)=Ψ-ψk/Ψ of the integrand function (a). Eigenvalues versus parameter ζ (b and c). Continuous line for exact values and dotted one for Improved Solutions: square tick for Linear (ILS) and bullet for Quadratic (IQS).
(a)
(b)
(c)Real Eigenvalues. It is expected that the real eigenvalues of (59) will be relatively close to the parameters μj in lightly damped systems, as described by Adhikari and Pascual [19]. Naming λ=λ(ζ) to the nonviscous eigenvalue associated with a particular parameter μr, then λ0=limζ→0⁡λ(ζ)=-μr. Indeed, multiplying (60) by λ+μr results in(63)λζ+μrλ2ζ+ωn2+2λζωnζNμr+2λζωnζN∑j=1j≠rNλζ+μrλζ+μjμj=0.Taking limits ζ→0 in the previous expression,(64)λ0+μrλ02+ωn2=0.Since λ0∈R for the nonviscous eigenvalues, λ02+ωn2>0 and consequently λ0+μr=0. Therefore, the nonviscous eigenvalue associated with μj can be expanded around the point λ0=λ(0)=-μj. 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 s=-μj (see (59)), the derivatives from (15) and (17) cannot be calculated. To avoid this problem, (63) with r=1 will be used. Thus, λ′(0) and λ′′(0) can be calculated taking derivatives with respect to ζ in the expression:(65)λζ+μ1λ2ζ+ωn2+2λζωnζNμ1+∑j=2Nλζ+μ1λζ+μjμj=0. Evaluating (65) at ζ=0, and after some operations, the following results are obtained: (66)λ′0=2ωnμ12ωn2+μ12,λ′′0=8ωn2μ13ωn2+μ122μ12-ωn2ωn2+μ12+∑j=2Nμjμj-μ1.The computation of the improved solutions from the differential equation (20) requires that G(s) is well defined at the initial point. Since λ0=-μj is a pole of G(s) and of ∂G(s)/∂s, the Taylor series up to the first order and second order will be used to estimate the eigenvalues.Figure6 shows results for a numerical example with ωn=3 rad/s and Biot’s model with N=3 kernels and μj={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,λIV0,…, from (65) can be extracted and new terms can be added.Figure 6
Exact solution and Taylor series approximation (quadratic and linear) of the nonviscous eigenvalues versus damping ratioζ. Biot’s model with N=3 kernels and nonviscous parameters μj={1,5,10} rad/s.


## 5.4. Multiple-Degree-of-Freedom Systems

In this last point, the application of the proposed method for multiple degrees-of-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 considern degrees-of-freedom (dof) vibrating system. The dofs’ time-domain response is represented by a vector u(t)∈Rn. With help of the Finite Element Method, the mass and the stiffness matrices of the system, denoted, respectively, by M,K∈Rn×n, 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 fdt∈Rn are proportional to the history of the dofs’ velocities via kernel hereditary functions. These functions are the entries of a matrix denoted by G(t)∈Rn×n, also assumed symmetric. Consider(67)fdt=∫-∞tGt-τu˙τdτ.Assuming small displacements, the motion equations linearly relate the dynamic balance among inertiaMu¨, elastic Ku, damping fd(t), and external forces fe(t). The general form of the system of integrodiferential equations can be written by (68)Mu¨+∫-∞tGt-τu˙τdτ+Ku=fet,u˙0=v0,u0=u0.Under these conditions, the modes of the system can be obtained as the nontrivial solutions of the free-motion problem obtained considering fe(t)=v0=u0=0 in (68). Thus, checking functions of the form u(t)=uest, we obtain (69)s2M+sGs+Ku≡Dsu=0,where G(s)=L{G(t)}∈Cn×n represents the Laplace transform of the viscoelastic kernels and D(s):C→Cn×n is thedynamic 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 (70)det⁡Ds=0.The complete set of eigenmodes of this problem present 2n+q distinct values distributed as(71)λ1,…,λn,λ1∗,…,λn∗,σ1,…,σq,where λj,λj∗∈C, 1 ≤j≤n, are n conjugate complex pairs, corresponding to the modes with oscillatory nature. The rest, σj∈R-, 1≤j≤q, are negative real numbers that represent overcritically damped modes callednonviscous 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 uj,uj∗∈Cn,1≤j≤n to the complex eigenvectors associated with λj,λj∗, and aj∈Rn,1≤j≤q to the eigenvector associated with the real eigenvalue σj.In the analysis and search of solutions of the eigenproblem, the proportionality ofn-dof system, that is, the modal decoupling capability, becomes of special importance. As known, the undamped problem obtained from G(t)≡0 can always be diagonalized in the modal space of matrices M and K. Denoting by ϕj∈Rn to the jth mass-homogenized real mode (in column) and by ωj to its associated natural frequency, the classical orthogonal relations can be written as (72)ϕjTMϕk=δjk,ϕjTKϕk=ωj2δjk,where δjk is the Kronecker delta. The modal matrix Φ groups in columns n undamped modes and is used to change the dof’s to the modal coordinates by u=Φq. Hence, (69) can be expressed as (73)s2IN+sΓs+Λq=0,where Λ=ΦTKΦ=diag⁡ωj2 is the diagonal squared natural frequency matrix and Γ(s)=ΦTG(s)Φ 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 calledproportional, since the necessary and sufficient conditions for the modal decoupling are directly related with proportional relationships between the damping matrix G(s) and the dynamic matrices M,K [31]. In some cases of nonproportionality, the matrix Γ(λj) is diagonally dominant but not purely diagonal, fact which is equivalent to assume as true: (74)∑l=1l≠knΓklλj<Γkkλj,∀1≤k,j≤n.The system is then said to belightly 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 (75)det⁡s2In+sΓs+Λ≈∏j=1ns2+sΓjjs+ωj2.Hence, the set of eigenvalues can be obtained from the following n decoupled equations: (76)Djs=s2+sΓjjs+ωj2=0,1≤j≤n,where Γjj(s)=ϕjTG(s)ϕj. Equation (76) represents the jth 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 Γjj(s,θ)=ϕjTG(s,θ)ϕj to refer to the jth modal damping function. For certain particular value of this parameter, say θ=θ0, (76) admits a closed-form analytic solution, denoted by λj(θ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 (77)λj1θ=λjθ0+λj′θ0θ-θ0,(78)λj2θ=λjθ0+λj′θ0θ-θ0+λj′′θ0θ-θ022,where the subscripts (·)j1 and (·)j2 denote the first- and second-order approximations of the jth eigenvalue, respectively. Expressions for λj′ and λj′′ 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 j: (79)Ψjs,θ=∂Γjj/∂θΓjjs,θ+s2+∂Γjj/∂sas a generalization of that of (21). It has been proved that theImproved Linear and Quadratic Solutions (ILS and IQS, see (24)) defined as (80)ILS:λ^j1θ=λjθ0exp⁡-∫θ0θΨjλj1ϑ,ϑdϑ,(81)IQS:λ^j2θ=λjθ0exp⁡-∫θ0θΨjλj2ϑ,ϑdϑ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 Figure7. 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 j=0,6 of fixed boundaries are related with zero dof’s so that u0=u˙0=u6=u˙6=0. Thus, the constitutive equations can be written as(82)Rj-1,j+TrαdαRj-1,jdtα=kΔuj+cTrαdαΔujdtα,j=1,3,4,5,where Δuj=uj-uj-1 and c,α, and Tr 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 c>1,0<α<1,Tr>0. The coefficient k, is the linear, static rigidity.Figure 7
Multiple degrees-of-freedom example: lumped-mass dynamical system with viscoelastic links based on the fractional derivatives.The kernel functionG(t) 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 G(s) can easily be calculated by simply applying to the fractional derivatives of (82) the Laplace transform and using its properties. If (83)Rj-1,js=LRj-1,jt,Δu^js=u^js-u^j-1s=LΔu^jtare the Laplace transform of the reactions and relative displacements, respectively, then (84)Rjs=k1+cTrsα1+TrsαΔu^js≡k+sGsΔu^js,where(85)Gs=ksc-1Trsα1+Trsα.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 Figure7, resulting in(86)s2M+sGs+Ku^s=0,where M=mI5 and (87)K=k2-1000-12-1000-12-1000-12-1000-12,Gs=Gs1000001-1000-12-1000-12-1000-11≡GsQ.The viscoelastic model is controlled by three parameters directly related to the damping, sayc, Tr, and α. In multiple-dof 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 Tr=0 and c=1 transform the nonlinear eigenvalue problem of (86) into the undamped linear one; that is, for these particular values, it is verified that G(s)≡0. Among these parameters we chose c 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 asG(s,c) (note that now our parameter is θ≡c). Solving the undamped eigenvalue problem we can extract the natural frequencies and the normal modes {±ωj,ϕj}j=15. Thus, the jth modal damping function is(88)Γjjs,c=ϕjTGs,cϕj=GsϕjTQϕj≡GsQj.After some straight operations the function Ψj(s,c) can be expressed as (89)Ψjs,c=∂Γjj/∂θΓjjs,θ+s2+∂Γjj/∂s=1+sTrααc-1+2s21+sTrα2/kQjsTrα.Moreover, considering the eigenvalues as c-dependent functions, the first- and second-order approximations based on the Taylor series expansion are (90)λj1c=λj1+λj′1c-1=iωj-c-1kQjTriωjTrα-121+iωjTrα,λj2c=λj1+λj′1c-1+λj′′1c-122=λj1c-ic-12k2Qj2iωjTr2α1-2α+iωjTrα8ωj31+iωjTrα3.As aforementioned, the current method improves these approximations proposing the so-named ILS and IQS solutions whose expressions for this particular case are (91)ILS:λ^j1c=iωjexp⁡-∫γ=1cΨjλj1γ,γdγ,(92)IQS:λ^j2c=iωjexp⁡-∫γ=1cΨjλj2γ,γdγ,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 jth eigenvalue can be written as sj=s0j+δ, with s0j being certain point adequately chosen and δ being an unknown value to be obtained. It is assumed that δ/s0j≪1; hence the solution is obtained expanding the characteristic modal equation Dj(s) from (76) around s0j up to the first or up to the second order. Thus, the first-order approximation yields (using the same notation as [18]) (93)sj1AP≈s0j+δ1=s0j-Djs0j∂Djs0j/∂s,while the second-order approximation yields(94)sj2AP≈s0j+δ2=s0j+-B±B2-4AC2A,where (95)A=12∂2Djs0j∂s2,B=∂Djs0j∂s,C=Djs0j.The lowest of the two roots (in absolute value) is chosen sincea priori the correct value is not known. Adhikari and Pascual proposed as initial guess s0j=-ξjωj+iωj1-ξj2, where ξj=lims→0⁡Γjj(s)/2ωj, 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 s0j=iωj, 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 thedistance 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 ofloss factor peak, ηm, which can be used as a measure of the level of damping [42–45]. According to the experimental evidence [44], values within the range 0.1<ηm<0.50 can be considered as moderate damping, while a loss factor 0.5<ηm<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.Table 2
Viscoelastic model based on fractional derivatives. Values of the damping parameters for the two considered damping cases.
Damping level
c
α
T
r (sec.)
η
m
Moderate damping (MD)
10
0,50
1
0
-
6
0,45
High damping (HD)
500
0,50
1
0
-
6
0,90The 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.Figure 8
Relative error in the computation of eigenvalues: proposed (solid lines) and Adhikari-Pascual [18] (dashed lines). (a) Moderate damping. (b) High damping.
(a)
(b)


## 6. 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 degree-of-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 thekth order Taylor series expansion presents O(hk+1), it is demonstrated that the current method presents O(hk+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 degrees-of-freedom systems.



---

*Source: 9634103-2016-01-06.xml*