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.
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, 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 bywhere 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 and are the initial position and velocity of the mass and is the external, time-dependent, applied force.
With reference to the mathematical form of , 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 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 . 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–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 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, 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.
2. Dynamics of Single Degree-of-Freedom Viscoelastic Systems
The characteristic equation of the eigenvalue problem can be obtained trying solutions with the form in the free-motion equation obtained from (2). ConsiderHere, as mentioned, is the Laplace transform of the kernel function . 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 toThe 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 eigenvalues, with . These can be arranged in a set with the formwhere are a pair of complex conjugate roots which induce oscillatory motion. The rest 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 formwhere 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.
3. 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 depends not only on the variable but also on a set of damping parameters which controls the viscoelastic behavior of the system.
In the general case, any eigenvalue 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 . Denoting as , let us assume now that there exists a particular value (or initial point) of , say , so that the characteristic equation can be solved analytically in terms of , obtaining the related initial eigenvalue : Under these conditions, the parameter 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 . The method is constructed assuming that the eigenvalue at , say , 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 if
Under this hypothesis, it is also guaranteed that there exist (i) an open interval with and (ii) an open neighborhood around the point in the complex domain , such that the function will be unique and continuously differentiable in and its graph verifies Notice that the damping function has been represented directly as in the implicit function theorem due to the fact that the parameters 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.
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 th order can be obtained around the initial point , provided that the derivatives evaluated at this point are available. ConsiderwithThe computation of 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 . 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 formwhere 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 . 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 , the solution will be named Improved Linear Solution (ILS). If the second-order approximation 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.
4. 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:(H1)If is the initial point verifying (8), then .(H2)Let and be the neighborhoods of and , 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 in the set . In addition, and its partial derivatives with respect to up to order 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 :
The hypothesis (H1) allows assuring the existence of the function in the interval . One of the theses of the implicit function theorem is the existence of the set 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 is a function from (21), then a real and positive number exists, , such that .
Property 2. Let be the eigenvalue as function of the parameter . Then for any integer there exists a real and positive number, , such that .
Notice that with the previous nomenclature . 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 . Assuming that hypotheses (H1), (H2), and (H3) are satisfied, there exist two positive real numbers, and , 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 .
First, using (H3) and the subtraction in Taylor series, the previous norm for iswhere 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), for ,Now the error can be calculated; using the definitions of both functions,
The theorem states that the approximation order of the improved solution is , while approximating with Taylor series up to the th term gives . 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, , where is the viscous damping coefficient. The Laplace transform of the kernel function becomes the -independent function . In general, instead of the damping ratio is commonly used. Both are related by the well-known expression , and consequently the mathematical model depends only on the single parameter . The mass-normalized characteristic equation (3) can be rewritten aswhere now . 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 or up to the quadratic one .
Taking as the initial point, the solution of the characteristic equation corresponds to that of the undamped system, . 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 () and IQS () can be computed solving the following differential equation:
Substituting in (38) the first order approximation of the root 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 from (37); the differential equation results in and the explicit expression of IQS can also be easily integrated, resulting in
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 . For this numerical example an undamped natural frequency 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 far from the initial value 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.
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 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 . Therefore, the nonviscous behavior will be controlled by , so that when the damping model is viscous with coefficient [34, 35]. Otherwise, when the behavior is highly nonviscous. The characteristic equation in the Laplace domain can be written as where . The coefficient of the limit viscous model can again be expressed in terms of the damping ratio . Under these conditions, the viscoelastic model can be written as function of the two parameters, and , as follows:
Introducing a new parameter with , the previous characteristic equation can be rewritten asEquations (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 . 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 byAs an initial point, any root of the equation can be taken. Without loss of generality we will use the positive imaginary part . From (15) and (17) the values of can be computed, obtaining after some operations the expressionsHence, 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 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 and have been plotted versus for an oscillator with natural frequency rad/s. Three damping ratios have been considered for comparison. In addition, the relative error of the exact integrand function ,is shown. Notice that the exact eigenvalue from (46) is required for the error computation.
(a)
(b)
(c)
Figure 3 shows that, as expected, close to the initial point , the proposed method eigenvalues are very accurate for all damping levels. Moreover, it can be seen that in the range , inequality that, as predicted from the theory of Section 4, implies that IQS is more precise than ILS. On the contrary, when , the errors are and, in addition, while increases monotonically, 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 (). The solution of this case is somewhat equivalent to the expansion of the eigenvalues in their Taylor series around or . 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 equationThe complex roots of this equation for 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 . 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 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 and reverse their relative position but now both tend to zero for high values of . Except close to the origin (), 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.
(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 aswhere are the relaxation parameters and 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 Note that the previous equation can be transformed into a polynomial of order . The set formed by eigenvalues can be separated, on one side, in a subset with the two complex conjugate eigenvalues and, on the other, in another with nonviscous eigenvalues that in general are real and 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 permits obtaining closed solutions for both the complex eigenvalues and the nonviscous ones.
Complex Eigenvalues. The value produces the eigenvalues of the undamped system in (60); that is, . Expanding the solution in its Taylor series around this point, the proposed method allows obtaining . The first and second derivatives, and , can be calculated from (15) and (17). After some operations, The second-order Taylor series become function of . 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 kernels. The relaxation parameters are taken as rad/s. Three different natural frequencies 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 and as before. The damping ratio has again been represented log-scaled. As in the previous case , (several orders of magnitude) , obtaining better approximations for IQS than those for ILS. For both viscous and nonviscous systems the agreement is very good, except for and rad/s (viscous behavior), where the slope of the imaginary part becomes almost vertical.
(a)
(b)
(c)
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 . Indeed, multiplying (60) by results inTaking limits in the previous expression,Since for the nonviscous eigenvalues, and consequently . Therefore, the nonviscous eigenvalue associated with can be expanded around the point . Without loss of generality, it can always be assumed that is associated with the first parameter . Thus, for any value , the Taylor series expansion is now . 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 will be used. Thus, and can be calculated taking derivatives with respect to in the expression: Evaluating (65) at , 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 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 rad/s and Biot’s model with kernels and rad/s. As expected, the eigenvalues show very good accuracy for damping range . Better results could be achieved through expansion of the Taylor series up to higher orders. For that, successive derivatives, , from (65) can be extracted and new terms can be added.
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 consider degrees-of-freedom (dof) vibrating system. The dofs’ time-domain response is represented by a vector . With help of the Finite Element Method, the mass and the stiffness matrices of the system, denoted, respectively, by , can be assembled. In general, are symmetric and positive definite and semidefinite, respectively. The viscoelastic damping is introduced in the system assuming that the dissipative forces are proportional to the history of the dofs’ velocities via kernel hereditary functions. These functions are the entries of a matrix denoted by , also assumed symmetric. Consider
Assuming small displacements, the motion equations linearly relate the dynamic balance among inertia , elastic , damping , and external forces . 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 in (68). Thus, checking functions of the form , we obtain where represents the Laplace transform of the viscoelastic kernels and 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 distinct values distributed aswhere , 1 , are conjugate complex pairs, corresponding to the modes with oscillatory nature. The rest, , , 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 to the complex eigenvectors associated with , and 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 can always be diagonalized in the modal space of matrices and . Denoting by 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 . Hence, (69) can be expressed as where is the diagonal squared natural frequency matrix and 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 and the dynamic matrices [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 Hence, the set of eigenvalues can be obtained from the following decoupled equations: where . 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 to refer to the th modal damping function. For certain particular value of this parameter, say , (76) admits a closed-form analytic solution, denoted by . 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 are where the subscripts and 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 of fixed boundaries are related with zero dof’s so that . Thus, the constitutive equations can be written aswhere 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 . The coefficient , is the linear, static rigidity.
The kernel function 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 inwhere and
The viscoelastic model is controlled by three parameters directly related to the damping, say , , 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 and transform the nonlinear eigenvalue problem of (86) into the undamped linear one; that is, for these particular values, it is verified that . 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 (note that now our parameter is ). Solving the undamped eigenvalue problem we can extract the natural frequencies and the normal modes . Thus, the th modal damping function isAfter 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 , with being certain point adequately chosen and being an unknown value to be obtained. It is assumed that ; hence the solution is obtained expanding the characteristic modal equation from (76) around 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 yieldswhere 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 , where , 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 , 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–45]. According to the experimental evidence [44], values within the range can be considered as moderate damping, while a loss factor 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.
(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 the th order Taylor series expansion presents , it is demonstrated that the current method presents .
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.
Appendices
A. Proof of Property 1
The function is well defined, provided that the denominator Hence, let be a positive real number such thatIn addition, according to hypothesis (H2), the function and its derivatives with respect to are bounded in the set . Thus, a positive real number exists such that Therefore, taking into account (A.2) and (A.3), the function can be bounded by the positive real number :
B. Proof of Property 2
The function is solution of the characteristic equation; therefore the following expression is verified: where the function . By hypothesis, the function is bounded in the set . Hence, there exists a real positive number such that , . Reordering (B.1) and taking absolute values, As a consequence, the function verifies the inequality The roots of the second-order polynomial are the real numbers . Therefore, the function is bounded by the interval because . Finally, the expression of the eigenvalue bound leads to
In order to bound the derivatives , , (20) relates the first with functions and whose bounds have already been calculated by (A.4) and (B.5). Thus, Higher derivatives of with respect to the parameter will be expressions function of the successive derivatives of the viscoelastic function which, from hypothesis (H2), are bounded in the interval . Consequently, once and are calculated, the numbers , for , can be obtained by induction.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.