Approximate method for temperature-dependent characteristics of structures with viscoelastic dampers

A method for analyzing the influence of temperature on the dynamic characteristics of structures with viscoelastic dampers is proposed in this paper. Dampers which are described by so-called fractional rheological models are considered. The temperature–frequency superposition principle is used to describe the influence of temperature on the dynamic characteristics. The concept of continuous dependence of damping on an artificially introduced parameter is adopted for solving, in an approximate way, the nonlinear eigenvalue problem from which dynamic characteristics are determined. The correctness and effectiveness of the method was verified by means of two examples.

steel frame model. The analysis was done in the time and the frequency domain, but the dimensionless damping ratios and natural frequencies were not determined.
It should be noted that changes in the temperature of the viscoelastic material in the damper are the result of external influences on the one hand and of the so-called self-heating phenomenon on the other. This phenomenon, consisting in the rise in temperature due to energy dissipation in the viscoelastic material, has been the subject of numerous studies. A number of coupled thermal-mechanical models have been developed, which are described in [8][9][10][11][12].
It is well known that temperature changes have a significant impact on the parameters of the damper model adopted [3][4][5][6]. The most commonly used models of viscoelastic dampers are classic rheological models described by the Prony series, fractional derivative models or power-law representations [10,13,14]. In recent times, fractional derivative models have often been used because of the small number of parameters needed to be identified. Model parameters are identified using the values of the so-called complex modulus, measured in experimental tests (for example, in dynamic mechanical analysis tests [14,15]). From the experimental results, the so-called master curves of the real part (storage modulus) and imaginary part (loss modulus) of the complex modulus are created. The values of the complex modulus for other temperatures are obtained by shifting the master curves [2,[14][15][16][17]. Shifting the master curves changes the model parameters, and this requires their re-identification. This paper presents an easy method for determination of the fractional derivative model parameters of the viscoelastic damper for different temperatures without their re-identification.
The main purpose of the paper is to develop a method for determining the basic dynamic characteristics, i.e., natural frequencies and dimensionless damping factors, for structures with viscoelastic dampers depending on temperature changes. The calculation of these characteristics requires the solution of a nonlinear eigenvalue problem. This type of eigenvalue problem can be solved using the mathematically advanced continuation method [18]. In [19], the method of approximating an eigenvalue problem solution was presented. It was extended in [20] to the case where the viscoelastic properties of the attenuators are described by fractional derivatives.
In order to determine the dynamic characteristics of a system with viscoelastic dampers for different temperatures, it is necessary to solve the eigenvalue problem repeatedly due to changes in the parameters of the damper models. However, on the basis of the method described in [20], a method of solving the nonlinear eigenvalue problem will be developed that takes into account the effect of temperature changes on the parameters of the damper models in a direct way. Details of the proposed method will be presented in this paper.
The remaining part of the paper is organized as follows. The fractional damper model is described in Sect. 2. The influence of temperature on the viscoelastic damper parameters is described in Sect. 3. The nonlinear eigenproblem is formulated in Sect. 4. The analytical solution of the nonlinear eigenproblem for proportionally damped systems is presented in Sect. 5. Section 6 describes the approximate method of solving the nonlinear eigenproblem. Section 7 presents an approximation of the dynamic characteristics of a structure according to the temperature changes in a direct way. An alternative method for approximating the dimensionless damping ratios of a structure is described in Sect. 8. In Sect. 9, numerical examples are presented. The final remarks are provided in Sect. 10.

Fractional damper model
Adopting a suitable viscoelastic damper model is a difficult problem: on the one hand, it should accurately describe the dynamic behavior of the damper; on the other, it is important to minimize the number of model parameters. The fulfillment of these two conditions can be facilitated by the use of fractional derivatives in the mathematical damper model.
The simplest fractional damper models are the fractional Maxwell model and the fractional Kelvin model. Each of them consists of the elastic element described by the stiffness parameter k, the Scott-Blair element described by the viscosity parameter c, and the parameter α, which expresses the order of the fractional derivative in the damper motion equation. These models are relatively simple from a mathematical point of view but in some instances they reproduce poorly the behavior of real viscoelastic materials. More complex models are a combination of the Maxwell and Kelvin models. In this paper, the model shown in Fig. 1 is taken into account. The considered model consists of the simple Maxwell and the simple Kelvin ones connected in parallel. Later in this paper, these two simple parts of the model will be called the Maxwell and Kelvin elements, respectively.
The total damper response force is the sum of the forces of the Kelvin and Maxwell elements, i.e.: where symbol u 0 (t) denotes the response force of the Kelvin element and u 1 (t) is the response force of the Maxwell element, respectively. The response force of the Kelvin element is described by the following formula: where D α t (·) denotes the fractional derivative of the α order with respect to time t. Here, the Riemann-Liouville definition of the fractional derivative is used. The symbol q (t) is the relative damper displacement: where q i (t) and q j (t) are the displacements of ith and jth nodes of the damper model (see Fig. 1). The response force of the Maxwell element is governed by the equation: where After using the Laplace transformation with zero initial conditions to Eqs. (2) and (4), the transform of total force of the fractional damper model response can be written as follows: whereū (s) and q (s) are the Laplace transforms of u (t) and q (t), respectively, and s = iλ is the Laplace variable. The symbol λ denotes the frequency of vibration. The constitutive equation of damper (5) can also be written using the complex modulus K (λ): A comparison of Eqs. (5) and (6) enables the storage modulus K (λ) and loss modulus K (λ) to be expressed using the damper parameters:

Influence of temperature on viscoelastic damper parameters
The complex modulus depends not only on the frequency but also on the temperature of a viscoelastic material. For most viscoelastic materials, the dependence of the complex modulus on temperature can be described by the so-called temperature-frequency correspondence principle [21]. If the storage modulus or loss modulus master curve, which is its plot in the frequency domain, is known for a certain temperature T 0 , called the reference temperature, then, according to this principle, it is possible to define the plot of this module for any other temperature T . This is done by shifting the graph at the reference temperature. According to the temperaturefrequency correspondence principle, the following relationship for the complex modulus can be written: where α T denotes the so-called shift factor and λ 0 is the reference frequency. In general, the horizontal and vertical shift factors must be used for the description of viscoelastic materials [2]. However, there are not enough experimental data concerning the vertical shift factor of the viscoelastic material used in dampers (compare [8,9,11,13,[22][23][24]). However, if the viscoelastic material can be classified as a thermorheologically simple material, then only the horizontal shift factor is required [23].
The value of the horizontal shift factor α T is most often determined using the empirical William-Landel-Ferry formula: where C 1 and C 2 are the empirical constants and T = T − T 0 . Moreover, it is assumed that the fractional free volume increases linearly with respect to time. An additional assumption is that, as the free volume of material increases, its viscosity rapidly decreases (compare [14]). The temperature-frequency correspondence principle for thermorheologically simple viscoelastic materials is illustrated in Fig. 2 and the following relationship can be written: The storage and loss moduli are described for the reference temperature T 0 by the formulae: For temperature T , the damper model parameters are different, and Eqs. (12) and (13) take the form: However, the moduli K (λ, T ) and K (λ, T ) for any temperature T can be derived from Eqs. (12) and (13) taking into account the temperature-frequency superposition principle. Replacing λ 0 by α T λ gives the relationships: The comparison of Eq. (14) with Eqs. (16) and (15) with Eq. (17) leads to the following results: whereα T = α α T .

Nonlinear eigenproblem
The equation of motion of structure with viscoelastic dampers could be written in the following form [18,20]: where M, C and K are the mass, structural damping and stiffness matrices of structure, respectively. The symbol q (t) denotes the vector of displacements in a time domain. Moreover, p (t) is the vector of external load and f (t) is the vector of interaction forces between the dampers and the structure (in short, the damping forces). The vector f (t) can be written in the form of the sum of vectors of forces induced by the individual dampers: where m is the total number of dampers. As a result of the Laplace transformation of Eq. (19) for the zero initial conditions, the following expression is obtained: Taking into account Eq. (5), the vector of damping forces can be written as: where L r is the location matrix of the r th damper. The equation of motion of structure with viscoelastic dampers in the frequency domain takes, finally, the following form: The basic dynamic characteristics, such as natural frequencies and dimensionless damping factors, are calculated after solving the appropriately defined eigenproblem. The nonlinear eigenproblem of the structure with viscoelastic dampers is obtained from Eq. (23) by subtracting the zero external load vector, i.e.: where Taking into account Eq. (18), it is possible to introduce the influence of temperature on the damper parameters into the equation of motion: where In a vast majority of cases, an analytical solution to the eigenproblem (26) is impossible. However, numerical methods have been developed (see, for example [18]). In this paper, the method of which the main assumptions were proposed by Lazaro [19] is presented and used.
The basic dynamic characteristics of structures are their natural frequencies ω k and dimensionless damping factors γ k , which are calculated from the formulae: where μ k and η k are the real and imaginary parts of the kth eigenvalue s k being the solution to the eigenproblem (26).

Solution to eigenproblem of proportionally damped system
In general, the nonlinear eigenproblem (28) has to be solved numerically. But in the special case of a structure in which the effect of dampers is proportional to the stiffness matrix, a "semi-analytical" solution can be given. The "semi-analytical" solution means that the characteristic equation is obtained but is nonlinear and has to be solved numerically.
If the system is proportionally damped, the eigenproblem (26) can be written in the following form: where κ denotes the proportionality coefficient. The coefficient κ fulfils the following relationship: After transformation, Eq. (31) can be written in the form of the linear eigenvalue problem: Taking into account thats = iω, the solution of the eigenproblem (31) comes to the solution of the equation: Above, the symbol ω denotes the natural frequency of the system without damping. The nonlinear characteristic Eq. (35) has to be solved numerically, for example, using the Newton method. In Eq. (35), the influence of temperature change on the damper model parameters can directly be taken into account by considering Eq. (18). Then, Eq. (35) takes the following form:

Approximate method of solving nonlinear eigenproblem
In the method proposed by Lazaro [19], the artificial parameter p is introduced to Eq. (26): where The parameter p is treated as a variable, and its value is in the range from 0 to 1. For p = 0, Eq. (37) describes the linear eigenproblem. For p = 1, the solution of nonlinear eigenproblem is obtained. It should be noted that s is a function of p. The first derivative of Eq. (37) with respect to p is as follows: where the symbol ∂D (s, p)/∂s denotes the derivative of D (s, p) with respect to p occurring as an explicit variable.
After premultiplication byq T ( p), Eq. (38) takes the form: The matrix D (s, p) is symmetrical. Thus, Eq. (37) can be written as: and from Eq. (40), the following expression is obtained: It is assumed that the eigenvectorq is scaled in the following manner: where a is a scalar. After taking into account the condition (43) and deriving the derivative of a matrix D with respect to p, Eq. (42) finally takes the form: where The next step of the presented method is to expand the D 1 (q, p) and D 2 (q, s, p) functions in the Taylor series in vicinity of p = 0: This operation leads to the following first-order ordinary differential equation: where Equation (49) is the nonlinear first-order differential equation which can be solved numerically. In this paper, the Euler method [25] is used. In order to calculate the magnitude dq dp p=0 , i.e., the sensitivity of the eigenvectorq to changes in the parameter p, Eq. (39) and the first derivative of Eq. (43) with respect to p are used. After some simple mathematical operations, the following matrix sensitivity equation is obtained: where

Analysis of temperature influence on dynamic characteristics of structures with viscoelastic dampers
Determining the dynamic characteristics of structures with viscoelastic dampers according to temperature changes within their specific range requires multiple solving of the nonlinear eigenproblem (26) because recalculation of the parameters of dampers for each temperature is required. This is quite troublesome because, as mentioned earlier, this task can only be solved by iterative methods. However, the use of the Lazaro assumptions [19] enables a direct determination of the course of the natural frequency and the dimensionless damping factor functions within a certain temperature range. Namely, the shift factorα T in Eq. (29) can be treated as a homotopy coefficient, as does the coefficient p in Eq. (38). After performing the corresponding mathematical operations described in Sect. 6, the following equation is obtained: where The next step is to expand the function D q, s,α T in the Taylor series forα T = 1: whereq 0 and s 0 are the eigenvector and the eigenvalue of the system, respectively, for the reference temperature T 0 α T = 1 . Finally, Eq. (61) takes the following form: where The sensitivities of the eigenvector and eigenvalue due to changes in the parameterα T are calculated from the following matrix sensitivity equation: where

Modal strain energy method
The modal strain energy method is used to approximate the dimensionless damping ratios of a structure [3,4]. The damping ratios are calculated from the following formula: where a k is the kth modal shape vector of the structure without damping. The symbols D I k and D Rk denote the imaginary and real parts of the dynamic stiffness matrix D (s) described by Eq. (27). In order to determine the matrices D I k and D Rk , the following expression has to be substituted into Eqs. (25) and (27): where ω k is the kth natural frequency of the structure without damping. After separating the real and imaginary parts of the dynamic stiffness matrix D (s), the matrices D I k and D Rk can be described as follows: The influence of the temperature changes on the parameter of the viscoelastic damper model can be taken into account using Eq. (18). Then, Eqs. (75) and (76) take the form: α T ω α c 0r cos απ 2 + k 1rαT ω α ν 1r sin απ 2 +α T ω α ν 2 1r + 2ν 1rαT ω α cos απ 2 + α T ω α 2 L r . (78) 9 Numerical examples 9.1 Example 1: the proportionally damped system In the first numerical example, the results obtained by the approximate method described in Sect. 6 are compared with the exact solution for proportional damping systems described in Sect. 5. The model of a four-story shear frame, i.e., the system with rigid beams and elastic columns, was used for the analysis. The system is shown in the schematic in Fig. 3a. The mass of the structure is lumped on the floor levels. The value of each lumped mass is m = 10,000 kg. The bending rigidity of the system is described by the story rigidity coefficients, k, which are also the same for each story (k = 1600 kN/m). Structural damping was omitted as irrelevant due to the nature of the study. In order to implement the damping proportionality condition, the dampers of the same parameters were located on all stories. The following parameters of the damper model were adopted: k 0 = 0.0 kN/m, c 0 = 0.0 kNs/m, k 1 = 10,000 kN/m, c 1 = 60 kNs α /m, α = 0.7. The nonlinear equation (35) was solved using the Newton method. In the calculation by the approximate method, it was assumed that the increment in the homotopy parameter, p, was p = 0.01. The nonlinear differential equation (49) was solved using the implicit Euler method. The results of the calculations are presented in Tables 1 and 2. For natural frequencies, the difference in the results is a maximum of 0.05%. In the case of dimensionless damping ratios, the maximum difference is 0.17%. It should be noted that the system with high damping is considered.
The example presented above shows that, although the method proposed by Lazaro [19] is an approximate one, its accuracy for systems with proportional damping, including systems with high damping, is very good at least in engineering applications.  Using the same model with proportional damping, the efficiency of the method of approximation of the dynamic characteristics of the system depending on the temperature changes, described in Sect. 7, was further analyzed. The temperature changes translate into changes in the coefficientα T . The range ofα T is from 0.4 to 2.0 in calculations, which corresponds to the range of α T between 0.3 and 2.7. The increment ofα T was 0.001. For example, for the material described in [10], the given range of α T corresponds to the temperature range of − 3.7-5.5 • C at the reference temperature T 0 = 0.2 • C. The results obtained from the numerical solution of Eq. (64) were compared with those calculated from Eq. (36). These comparisons are illustrated in Figs. 4 and 5. For proportionally damped systems, convergence of the results of the approximate method with those obtained from the exact solution can be considered quite good: the maximum difference in results for the natural frequencies is 0.63%, and that for the dimensionless damping factors is 3.81%.
In addition, the plots of dimensionless damping ratios calculated using the modal strain energy method are presented in Fig. 5. It can be seen that the convergence of the results obtained with this method with the exact solution is worse than that for the method proposed in Sect. 7 for the proportionally damped system. 9.2 Example 2: the non-proportionally damped system The effectiveness of the approximation method of calculating the dynamic characteristics of the nonproportionally damped systems according to temperature changes was verified by performing numerical calculations for the shear frame model of the building structure with ten degrees of freedom and six dampers. The system is shown in the schematic in Fig. 3b. It was assumed that changes in the temperature of the viscoelastic material affect the changes shift factor, α T . Calculations were performed twice for the same α T interval. In the first series of calculations, the damper model parameters were changed for the subsequentα T values according to Eq. 49, and then the problem was solved by the method described in Sect. 6. In the second series of calculations, the method described in Sect. 7 was used, where the subsequentα T values were substituted for Eq. (64). As in Example 1, the parameterα T was between 0.4 and 2.0, which translates into α T of 0.3-2.7. The limits of theα T range were determined in such a way as to obtain a satisfactory consistency of the results. The increment ofα T for the first calculation series was assumed to be 0.01. For example, for the material described in [10], the given range of α T corresponds to the temperature range of − 3.7-5.5 • C at the reference temperature T 0 = 0.2 • C. For the second series, the value 0.001 was used to increase the accuracy of the calculation. It should be emphasized that despite the decrease in the increment value ofα T , the speed of computer-aided calculations of the second series was much higher. For the first three natural frequencies, the maximum differences between the results were 0.45, 1.31 and 2.54%, respectively. For the dimensionless damping factors of the first three vibration modes, the effects were a bit worse. Namely, the maximum differences in the results were 7.69, 7.87 and 9.36%, respectively. The natural frequency functions versus α T are shown in Fig. 6. The graphs of the dimensionless damping ratios are given in Fig. 7. The values obtained from the modal strain energy method are also shown in this figure. It can be seen that the results obtained using all the presented methods are very similar.
Comparing the results for proportionally and non-proportionally damped frames, it is obvious that results for the proportionally damped frame are noticeably more accurate.

Concluding remarks
Viscoelastic dampers are often used nowadays to reduce the excessive vibrations of structures. The modal analysis in the frequency domain needs an eigensolution to the nonlinear eigenvalue problem resulting from the equation of motion for a damped, free vibration problem. Although exact solutions to the above-mentioned eigenvalue problem can be obtained numerically, approximated expressions for frequency of vibrations and non-dimensional damping ratios are very valuable in reducing the computational effort. It is more important in the case of systems with viscoelastic dampers because their parameters depend on temperature.
In this paper, the method of analysis of the influence of temperature on the dynamic characteristics of structures with viscoelastic damper is presented. The method enables the calculation of the sequences of natural frequencies and dimensionless damping ratios of a structure in a certain temperature range without the need of multiply solving the nonlinear eigenvalue problem. The obtained dynamic characteristics values are approximate. The degree of consistency of the results with those obtained by directly solving the eigenvalue problem for each set of damper's parameters depends on the considered range of temperature.
Significant increase in the efficiency of the proposed approximation method in relation to the conventional eigenvalue problem solution is due to the fact that instead of multiple complex matrix operations, the ordinary differential equation (64) is solved only once. In addition, the effect of temperature changes is directly taken into account when solving this equation. In the conventional approach, the eigenvalue problem would have to be solved many times-for each temperature.
The effectiveness of the presented method was verified numerically for two examples of structures with four and ten degrees of freedom, respectively. Several classic and fractional damper models were used to   Plots of three first natural frequencies for Example 2: case 1-the values calculated using the presented method taking into account the temperature changes directly, case 2-the values calculated for the case when the parameters of the damper model were changed each time before solving the eigenproblem Fig. 7 Plots of dimensionless damping factors of three first vibration modes for Example 2: case 1-the values calculated using the presented method taking into account the temperature changes directly, case 2-the values calculated for the case when the parameters of the damper model were changed each time before solving the eigenproblem, case 3-the values calculated using the modal strain energy method describe the properties of viscoelastic dampers mounted on structures. The results of computation show that the suggested method can be successfully applied to proportionally and non-proportionally damped structures. The relative error of approximately calculated eigenvalues depends on the damping level and on the range of considered changes of temperature. The results obtained using the proposed method were also compared with those obtained by the modal strain method, and the proposed method was found to lead to smaller relative errors.