Abstract

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

1. Introduction

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

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

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

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

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

2. Viscoelastic Constitutive Relationships

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

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

3. Equations of Motion in the Frequency Domain

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

4. Proposal of the Viscous Model

4.1. Complex Eigenvalues

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

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

4.2. Single-Degree-of-Freedom Systems

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

Figure 1 

The original single-degree-of-freedom model with viscoelastic behavior based on the fractional derivatives and its equivalent linear viscous model. denotes the -order fractional derivative.

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

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

4.3. Multiple-Degrees-of-Freedom Systems

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