Analysis of dynamic characteristics of viscoelastic frame structures

Planar frame structures made of a viscoelastic material are considered in the paper. The technically very important structures made of a homogenous material are contemplated. A family of rheological models (classic and fractional) are used to describe the mechanical properties of the viscoelastic material. In particular, the dynamic characteristics of the structures are of interest. A numerically very efficient method is proposed to determine such characteristics. The method requires the solution to the linear eigenvalue problem for corresponding elastic structures and the solution to a nonlinear, algebraic equation. The presented method is much more efficient than other methods where, very often, the continuation method is used to solve the nonlinear eigenvalue problem. The influence of temperature changes on dynamic characteristics is analyzed using the frequency–temperature principle. The results of several parametric analyses are presented and discussed. For the first time, the influence of temperature on the dynamic characteristics of beams has been studied in detail.


Introduction
Viscoelastic (VE) materials are often used to build dampers, layers of structural elements and whole structures in order to increase the ability of systems to dissipate more energy, in comparison with conventional ones [1,2]. The mechanical properties of VE materials depend mainly on the frequency of excitations and temperature [2][3][4]. A family of rheological models which take into account the dependence of the mechanical properties of such materials on excitation frequencies are proposed and can be found in several papers and monographs [3,[5][6][7][8][9]. Both the classic rheological models such as the Zener model and the generalized Maxwell model [5,7] and the so-called fractional models [6,8,9], in which the fractional derivatives are utilized, have been proposed in the past. The influence of temperature changes on the properties of VE materials is relatively rarely taken into account [10][11][12][13]. In those papers, empirical relationships derived on the basis of experiments were used [10] or the so-called frequency-temperature superposition principle was utilized [12,13].
Natural frequencies, non-dimensional damping ratios, modes of vibration and frequency response functions are very important characteristics of each mechanical structure. For elastic systems, these properties can be determined using well-known methods, described in many monographs (see, for example, [5,14]).
Methods which are devoted to determination of the dynamic properties of VE structures or elastic structures with VE dampers also exist (see, for example, [15][16][17][18][19][20]). For systems with VE elements modeled with the help of classic rheological models, the dynamic characteristics are obtained after solving the quadratic eigenvalue problem [21,22] which can be transformed to the linear eigenvalue problem using the state variables, as shown, for example, in [5,21,23]. The dimension of the above-mentioned linear eigenvalue problem is much greater than the quadratic one and, many times, it is greater when the state variable concept is used in the analysis of systems with elements described using the fractional derivatives (see [24]). Moreover, the physical insight into the problem can be lost due to the introduction of the state variables. Singh [20] proposed a method based on the Taylor series expansion of transcendental matrices appearing in the nonlinear eigenvalue problem and combined with Newton's eigenvalue iteration method. Asymptotic numerical techniques are used in [25,26] to determine the natural frequencies and the loss factors of VE sandwich structures. Also, Adhikari and Pascual [19,27] proposed iterative methods based on the expansion in the Taylor series of a part of the eigenvalue problem corresponding to the VE elements of the dynamic system at hand. A recursive procedure based on fixed-point iteration which always converges to the complex eigenvalues was proposed by Lazaro et al. [28]. The iterative method is proposed in [29] to solve the nonlinear eigenvalue problem for a non-viscously damped structure. The continuation method is proposed to solve the nonlinear eigenvalue problems for systems with VE dampers (in [16]) and VE layers (in [15]) of which the properties are described with the help of fractional derivatives. The homotopy method, a version of the continuation method, is applied in [13] to determine the natural frequencies and the damping ratios for the planar frame with VE dampers, modeled by the fractional rheological models.
The main aim of this paper is to present a very simple method for determination of natural frequencies, damping ratios and the frequency response functions (FRF). The method is applicable to a restricted, though technically very important, class of structures, i.e., to framework structures built of a homogenous VE material. The material can be described with a freely chosen rheological model (classic or fractional), and the parameters of that model must be identical for all the elements of the structure. In particular, a new method for solving the nonlinear eigenvalue problem, of which the computational efficiency is very high, is presented. In the method, it is necessary to know the solution to the linear eigenvalue problem for some elastic structures and the solution to a single algebraic equation to determine the dynamic characteristics of the VE structures.
The paper is organized as follows. In Sect. 2, the equation of motion of VE rods treated as the continuous system is derived. The characteristic equation derived in the Laplace domain is solved. Then, in Sect. 3, two finite element formulations for VE structures are presented. In the first formulation, a typical two-node element with the Hermite's shape functions is shown and the dynamic version of the finite element is described in the second formulation. In the first formulation, the problem of determination of dynamic characteristics is reduced to the solution to the linear eigenvalue problem and a single nonlinear algebraic equation for each eigenvalue. In the second one, the solution to the nonlinear eigenvalue problem for the corresponding elastic structure and the solution to a set of single nonlinear algebraic equations is needed. In both cases, the dimension of the eigenvalue problems is as for the elastic structures. Some remarks about the characteristic features of the proposed method and the method of calculation of FRF are presented in Sect. 4. How the frequency-temperature superposition principle could be applied to take into account the influence of temperature on the dynamic characteristics of considered structures is shown in Sect. 5. The results of several parametric analyses are presented and discussed in detail in Sect. 6. Finally, concluding remarks are formulated in Sect. 8. Moreover, several useful formulae are given in Appendices.

Equations of motion and their solutions for VE beams treated as continuum systems
The Euler theory is used to describe the dynamic behavior of beams made of VE materials. Some of the well-known equilibrium equations of beams are: where M(x, t), T (x, t) and N (x, t) are the bending moments, the shear and normal forces, respectively; u(x, t) and w(x, t) are the axial and transversal displacements of neutral axis of beams; p x (x, t) and p y (x, t) are the axial and transversal excitation forces per unit length of a beam. The dot indicates differentiation with respect to time t and (·) ,x = ∂(·)/∂ x.
The axial strain of the neutral axis of a beam ε(x, t), the angle of rotation of the beam cross section ϕ(x, t) and the beam curvature κ(x, t) are defined as: The constitutive relationship for the VE material written at a point level is in the following form: where σ x (x, z, t) and ε x (x, z, t) are the normal stress and the axial strain in the material fiber of which the coordinates are (x, z). G(t − τ ) is the kernel function, also known as the heredity function. From the mathematical point of view, any expressions which make the dissipation energy nonnegative function could be used as the kernel functions G(t − τ ) [30]. At the level of the beam's cross section, the constitutive Eq. (5) could be rewritten in the following way: where A and J are the area of the beam's cross section and the second moment of the cross section, respectively. The damping kernel functions of the non-viscously damped systems are usually defined in the frequency domain. Therefore, Eq. (6) can be equivalently expressed in the Laplace domain with zero initial conditions: where s is the Laplace coordinate andḠ(s) is the Laplace transform of the heredity function. The particular form of the sḠ(s) function for the most popular rheological models of VE material is given in "Appendix A".
The Laplace transformation of Eqs (1)-(4) results in: After using the non-dimensional variable ξ = x/l (l is the length of the beam span) and substituting relationships (7.2) and (10.3) into Eq. (8.1), the following equation of axial vibration in the Laplace domain is obtained: where Proceeding in a similar way, i.e., introducing relationships (7.3), (9) and (10.2) into (8.2), the equation of transversal vibration in the Laplace domain is obtained: where It should be noted that the form of Eqs (11) and (13) is formally identical with the well-known equation for the elastic beam. It is a direct consequence of the elastic-viscoelastic correspondence or the elastic-viscoelastic analogy [4]. This fact was also noted for the cantilevered VE (Kelvin model) beam in [31]. The solutions to the elastic beam are well known and described in detail in many monographs. In this paper, for the reader's convenience, only the most important relationships are quoted. Solutions to Eqs (11) and (13) are in the following forms: w(ξ ) =D 1 sin λξ +D 2 cos λξ +D 3 sinh λξ +D 4 cosh λξ, (16) respectively. The constants D 1 , D 2 ,D 1 ,D 2 ,D 3 andD 4 are determined from the appropriate beam's boundary conditions. The conditions of existence of non-trivial solutions result in the characteristic equations with respect to λ which take, in the case of axial vibration, the following form: In the case of transversal vibration and, for example, the fixed-fixed beam, the characteristic equation is: The solution to the characteristic equation results in a set of characteristic values λ k where k = 1, 2, . . . , ∞. Having the characteristic values λ k , it is possible to determine the corresponding set of characteristic Laplace numbers s k from Eqs. (12) or (14) for both considered cases (i.e., the axial or the transverse vibration, respectively). It constitutes the nonlinear algebraic equations of which the particular form depends on the rheological model of the beam's material. Taking into account Eq. (A.18), these equations are rewritten below in the following form: for the axial and transverse vibrations, respectively, and where ω e is the corresponding natural (axial or transversal) vibration frequency of the elastic beam under consideration with the Young's modulus E 0 . However, it is well known that in the case of axial and transverse vibrations, respectively. Therefore, relationships (19) and (20) could be written in a uniform way as the following equation: The number of solutions of Eq. (22) depends on the model of the beam's VE material. In the case of the classic and fractional Kelvin model, two complex and conjugate numbers are obtained. For the classic Zener model, Eq. (22) has three solutions, i.e., two complex and conjugate numbers and one real solution. However, it is found in many cases that only two complex and conjugate solutions exist for the fractional Zener model. The nonexistence of real eigenvalues for the fractional Zener model, used to describe the VE material of a multilayer beam, was proved by Lewandowski and Baum [15]. The same problem was discussed in paper [16] for the fractional Maxwell model. In the case of the generalized Maxwell model, we obtain (2 + m) solutions to Eq. (22), i.e., two complex and conjugate numbers and m real numbers. In order to fulfill the stability criteria of motion, the real solution and the real part of complex solutions must be negative numbers. The complex characteristic numbers s are connected with the oscillating solution to the equation of motion, while the real numbers s are connected with the non-oscillatory solution to the equation of motion. An interesting conclusion from the presented analysis is that there are infinitely many real solutions, provided they exist. This connection between the complex and real eigenvalues is found for the first time. Moreover, if the system under consideration has some repeated complex eigenvalues, then they are also repeated real eigenvalues. Moreover, it is easy to recognize that all eigenvectors (the mode of vibrations) are the real vectors for all considered rheological models of beam's VE materials. These modes of vibration are identical with those of the elastic beam of which Young's modulus is E 0 . All material points of the beam simultaneously reach their maximum and zero displacement, respectively. It is also important to underline that the consequence of such properties of the modes of vibration is the validity of the orthogonality conditions, which is well known for an elastic beam. This also means that normal coordinates can be used in a usual way to derive the uncoupled equations of motion. Having the complex solution s = μ + iη to Eq. (22), it is possible to calculate the natural frequency ω and the non-dimensional damping ratio γ of the VE beam using the following formulae: Two methods were used to solve Eq. (22). We used the procedure solve from MAPLE and the Newton method described briefly below. When the Newton method is used, a good initial approximation for s is required because of the local convergence of the Newton method. In this paper, s (0) k = ±iω ek is chosen as the initial approximation of s k where we try to determine the complex solution to Eq. (22). (i = √ −1 is the image unit, and ω ek is the k-th natural frequency of the elastic beam.) Having the i−th approximation s (i) of the complex solution to Eq. (22), the next approximation is obtained from the following formula: where f ,s (s) = 2s + ω 2 e θ ,s (s) is the first derivative of f (s) with respect to s. The derivatives of θ(s) for the considered rheological models are given in "Appendix A". The iteration process is finished when the appropriate convergence criteria are fulfilled.
This procedure must be somewhat modified if the real solution to Eq. (22) is needed. The real solution exists only for a VE material modeled with the help of the classic Zener model or the generalized Maxwell model. Since the classic Zener model can be understood as the generalized Maxwell model with one Maxwell element, only the generalized Maxwell model will be considered below. It is known that the real solution is in the vicinity of s = −1/τ . Assuming that the j-th real solution must be calculated, Eq. (22) is first multiplied by (1 + τ j s) and rewritten as follows: Now the Newton method is used, as usually, starting from initial approximation s (0) = −1/τ j . The corrected approximation is calculated using the following formulae: where and the iteration process is finished when the appropriate convergence criteria are fulfilled. For the classic Zener model, the function g(s) and its derivative are: In some cases, it is useful to know how the dynamic characteristics of beams depend on some design parameters p; for example, the E ∞ /E 0 or τ could be the above-mentioned parameters. In that case, the s( p), ω( p) and γ (p) curves can be determined using the incremental-iterative procedure, where it is assumed that the solution to Eq. (22) is known for a particular value of p = p i . Next, the value of the design parameter is incremented by p (i.e., p i+1 = p i + p) and the solution to Eq. (22), i.e., s( p i+1 ) is determined using the Newton method. After that, it is easy to find ω( p i+1 ) and γ (p i+1 ). In this way, the so-called ordinary points on the above-mentioned curves can be obtained. However, as will be demonstrated later, the bifurcation points may also exist on these curves for some VE models. At the bifurcation point, of which the coordinates are denoted as s b and p b , the two conditions written below must be satisfied: The coordinates of the searched bifurcation point are obtained from the two equations above. In particular, for the classic Zener material and when p = E ∞ /E 0 , Eqs. (30) and (31) can be rewritten in the form: which can easily be transformed to the following equations: Only real solutions of Eq. (33a), providing a positive ratio of the E ∞b /E 0 modulus, are of interest.
To end with the description of the Newton method, some features of the function f (s) must be discussed. It is complex argument function which is non-holomorphic in a general case. However, it is a holomorphic function if the classic Kelvin model is used in the description of VE material. For the classic Zener model and the generalized Maxwell model, the function f (s) is non-holomorphic but Eq. (22) can be transformed to the form where the holomorphic function will be on the right-hand side of Eq. (22). That version of Eq. (22) is obtained after multiplying Eq. (22) by (1 + s b τ ) (in the case of the Zener model) or by m i=1 (1 + sτ i ) (in the case of the generalized Maxwell model). For the fractional rheological models, f (s) is always a nonholomorphic function. Strictly speaking, for the non-holomorphic function, the first derivative with respect to complex argument does not exist (c.f. [32]) and the Newton method cannot be used to solve the equation f (s) = 0. In spite of these, in this work, Eq. (22) is solved with success with the help of the Newton method and the derivative f ,s (s) was calculated as usual, i.e., in the case of the real argument function. This approach, though used to solve the nonlinear eigenvalue problem, was also utilized with success in [15,16]. However, it was also discovered that the Newton method, as presented above, failed (for α = 1) in the vicinity of the bifurcation point where the considered mode of vibration becomes overdamped.

Formulation using the ordinary finite element
The finite element method is a general and widely accepted method used to solve many dynamic problems. Here, the planar, two-node rod element with six nodal parameters and the Hermite's shape functions and shown in Fig. 1 is used to describe the dynamic behavior of VE beams.
The vector of nodal parameters (in a local coordinate axis) is q where u 1 , u 2 , w 1 , w 2 and ϕ 1 , ϕ 2 are the axial nodal displacements, the transverse nodal displacements and the angles of rotation of the nodal cross sections, respectively. The vector of displacements w e (x, t) = [u(x, t), w(x, t)] T is written in terms of nodal parameters in the following way: where is the matrix of shape functions. The symbols , H 6 (x) denote the linear and the cubic Hermite shape functions, respectively. The vector of generalized strains e(x, t) = [ε(x, t), κ(x, t)] T depends on the nodal parameters as follows: where The vector of excitation forces is

and the vector of the inertia forces
After introducing the vector of generalized stresses s( T , the matrix form of the constitutive relationships (6) is: where in (42) In the Laplace domain, Eq. (39) takes the following form: where s is the Laplace variable, the quantities with the bar denote the Laplace transform of the corresponding quantities without the bars (for example,s(x, s) is the Laplace' transform of s(x, t)) and Using the above notations, the virtual work equation could be written as follows: where the symbols δw e (x, t), δq e (t) and δe(x, t) denote the vector of virtual displacements, the vector of virtual nodal displacements and the vector of virtual generalized strains, respectively. Moreover, T is the vector of nodal reactions and l e is the length of the finite element. In the Laplace domain, the virtual work equation takes the following form: Proceeding in the usual way, the second and third terms of Eq. (44) can be rearranged in the following way: whereP e is the Laplace transform of nodal excitation forces and M e is the well-known elemental mass matrix of the rod. Moreover, the right-hand side of Eq. (44) could be written, using the relationships (39) and (44) as follows: where the elemental VE stiffness matrix is It is useful to divide the above VE stiffness matrix into two matrixes: the elemental elastic stiffness matrix K e and the second matrix K ev which will be called the elemental viscous matrix. Precisely speaking, the last matrix expresses the damping properties of the rod only when the classic rheological models are used in the constitutive Eqs. (39) and (42). Taking into account relationship (A.18), which is valid for all of the considered rheological models, it is possible to writẽ where is the well-known elemental stiffness matrix of the elastic rod made of an elastic material of which the elastic module is E 0 . Moreover, E 0 = E 0 A 0 0 J and it is important to note that both of the matrices mentioned above are proportional and the factor of proportionality is θ(s). Finally, from the virtual work Eq. (44), the following equation for the finite element is obtained. Applying the usual finite element procedure, the equation of motion in the Laplace domain can be written in the following form: whereq andP are the global vectors of the Laplace transform of nodal displacements and excitation forces, respectively. Moreover, M and K are the global mass matrix and the global stiffness matrix of the elastic beam, respectively. In this paper, it is assumed that the whole structure is made of a homogenous VE material and one rheological model is used to describe the mechanical properties of that material. The dynamic properties of the VE structure are obtained assuming, firstly, thatP = 0, and solving the nonlinear eigenvalue problem where s is the eigenvalue andq is the eigenvector.
In a general case and when the properties of VE elements of the structure are modeled with the help of classic rheological models, the modes of vibration are complex and are also called the elastic modes of vibration, except that, for some rheological models (for example, the Zener and the generalized Maxwell model), the real modes of vibration, also called the non-viscous modes, exist. However, the real modes of vibration do not exist for the fractional Zener model [16].
It is difficult to solve the nonlinear eigenvalue problem of this type. Fortunately, in the considered case, the eigenvalue problem (25) can easily be transformed into the following linear eigenvalue problem whereq l is the eigenvector and ω 2 e the eigenvalue of the linear eigenvalue problem (53), then linked with the eigenvalue of the nonlinear problem by where ω e is the natural frequency of the elastic rod made of a material with the elastic modulus E 0 , and then solving the following equation with respect to s. Equation (54) is identical with Eq. (22) derived for the continuous model of the rod. Many well-known procedures, described, for example, in [5,14] can be used to solve the eigenvalue problem (53).

Formulation using the dynamic finite element
The dynamic analysis of rods can also be formulated using the so-called dynamic finite element method. As shown in Sect. 2, the modes of vibration of the elastic and VE rods are identical. In the dynamic finite element method, these modes are used as the shape functions. A brief description of the dynamic finite element formulation for the VE rods will be presented below. As in the previous subsection, the typical finite element is considered. The axial and transverse vibrations are analyzed separately. In both cases, the non-dimensional spatial variable ξ = x/l e is utilized.
If the axial vibration of the VE rod is considered, the mode of vibrationū(ξ ), their first derivative with respect to the spatial variable ξ and the Laplace transform of the axial force are: where D 1 and D 2 are constant and At the ends of the finite element, the following boundary conditions must be fulfilled:ū(0) =ū 1 ,N (0) = −R 1 ,ū(1) =ū 2 andN (1) =R 2 . These conditions lead to the following relationships: Now, it is easy to recognize that the stiffness matrix for the rods executing axial vibration is: where F 7 (λ e ) = (l 2 e /i 2 e )λ e /tgλ e , F 8 (λ e ) = λ e l 2 e /(i 2 e sin λ e ). When the VE elements execute transverse vibration, the mode of vibrationw(ξ ), their first derivative with respect to ξ , the Laplace transform of bending moments and shear force are: The first set of boundary conditions which must be fulfilled is:w(0) =w 1 ,w ,ξ (0) = l eφ1 ,w(1) =w 2 , w ,ξ (1) = λ eφ2 , and leads to the following matrix equation: where S = sin λ e ,S = sinh λ e , C = cos λ e ,C = cosh λ e . The second set of boundary conditions which must be satisfied is: These conditions provide another matrix equation in the form: From (67), we obtain d e = A −1 eq e which, when introduced into (69), leads to the relationshipR e = B e A −1 eq e . It means that the searched dynamic stiffness matrix is The inverse matrix to A e is in the following form: Taking into account relationships (70), (72) and (A.18), from (71), we obtain the dynamic stiffness matrix in the following explicit form: where Again, it is visible that the dynamic stiffness matrixK e is a sum of two parts: the stiffness matrix of an elastic rod and the VE stiffness matrix.
For the readers' convenience, the dynamic stiffness matrix of the rod executing both the axial and transverse vibrations is given in "Appendix B".
Applying the finite element procedure, the following nonlinear eigenvalue problem can be derived: whereq are the global vectors of the Laplace transform of nodal displacements and K(λ) is the global dynamic stiffness matrix of the elastic frame structure corresponding to the VE structure under consideration. As in the previous case, it is assumed that the whole structure is made of a homogenous VE material and one rheological model is used to describe the mechanical properties of that material. The symbol λ is the eigenvalue of the nonlinear eigenvalue problem for the above-mentioned elastic structure and is linked with the eigenvalue s of the VE structure through relationship (66) which is now reinterpreted and rewritten in the following form: wherem r ,l r andJ r are referred to the reference finite element. Additional details are given in "Appendix B". A non-trivial solution to the nonlinear eigenvalue problem (75) exists if det K(λ) = 0.
This condition is the basis of the Wittrick-Williams procedure for solving the eigenvalue problem (75) as presented in [33].

Remarks on modes of vibration, the orthogonality conditions and dynamic characteristics of frame structures
In a general case, when the properties of the VE elements of a structure are modeled with the help of classic rheological models, the modes of vibration are complex and called also the elastic modes of vibration [34]. The exception is that, for some rheological models (for example, the Zener and the generalized Maxwell model), the real mode of vibration (called also the non-viscous mode) exists.
It is easy to recognize that the linear eigenvalue problem (53) is identical with the eigenvalue problem from which the natural frequency of the considered structure is calculated but is made of an elastic material with the elastic modulus E 0 .
The matrices K and M in Eq. (53) are positive definite ones, which means that the eigenvalue ω 2 e is a real, positive number and the eigenvectorq l is a real vector. Moreover, sinceq =q l , the eigenvector of the nonlinear eigenvalue problem (52) is also a real vector. In addition, two eigenvectors of this kind (say,q i and q j ) fulfill the well-known orthogonality conditions: where δ i j is the Kronecker symbol. However, the eigenvalues of the above-mentioned problems (52) and (53) are different because s must be determined from Eq. (54) which is a nonlinear algebraic equation and, for a given ω e , has more than one solution. A number of solutions depend on the rheological model used to describe the beam's material. For example, for the Zener model, numerical experiments show that Eq. (54) has two complex, conjugate solutions in the case α = 1 and three solutions (two complex, conjugate solutions and one real solution) when α = 1.
In the case of elastic structures with n-th degree of freedom, we have n eigenvalues and eigenvectors but there are more than n solutions to the eigenvalue problem for VE structures because of the fact that, for one linear eigenvalue, more than one solution to Eq. (54) is obtained. Moreover, according to the stability conditions of motion, negative real eigenvalues and negative values of the real part of complex eigenvalues are required.
For example, when the fractional Zener model (α = 1) is used, we have 2n eigenvalues but in the case of the classic Zener model (α = 1), we have 3n eigenvalues (2n complex, conjugate eigenvalues and n real eigenvalues). Because the number of eigenvectorsq andq l is n, we conclude the vectorq is in fact a repeated eigenvector corresponding to a few (two or three in the case of the Zener model) eigenvalues s, determined from Eq. (54) for a given value of ω e . One consequence of this, for the considered type of VE structures, is that the orthogonality conditions do not hold for two eigenpairs of which the eigenvalues fulfill Eq. (54) for a chosen natural frequency ω e of the corresponding elastic structure.
For the VE beams and frame structures considered in this case, one above-mentioned mode of vibrationq l can be understood as the viscous mode (when it is an eigenpair with the complex eigenvalue s) or as the nonviscous mode (when it is an eigenpair with the real eigenvalue s). Moreover, it must be remembered that the proof of the orthogonality conditions of the eigenvalue problem (53) is derived using eigenvalues λ rather than s.
Having a complex solution s = μ + iη to Eq. (54), it is possible to calculate the natural frequency ω and the non-dimensional damping ratio γ of the VE structure using the following formulae: Often the frequency response functions (FRF) are used to characterize the dynamic behavior of structures. The FRF matrix Z(λ) is obtained assuming s = iλ and inverting the matrix D(λ) = [1 + θ(iλ)]K − λ 2 M. For the considered type of structures, this can be done either in a direct way or using a concept similar to that of normal coordinates as described below. For large systems, the second approach is more efficient.
The solution to Eq. (51) for s = iλ is assumed in the form: where a j is the eigenvector fulfilling Eq. (53) andx j (λ) is the Fourier transform of the normal coordinate x j (t) of the elastic structure. After introducing relationship (80) into Eq. (51) multiplied by a T i and taking into account the orthogonality conditions (78), the following relationship is obtained: The searched FRF matrix can now be written in the form:

Influence of temperature on dynamic characteristics of structures
The properties of VE materials significantly changed with temperature [4]. Temperature in VE materials can vary due to changes of the environmental temperature or due to the self-heating process in which the energy of vibration is dissipated and transformed into heat. The self-heating process is theoretically analyzed in [35] where it was found that, after a sufficiently long time, temperature in the solid VE damper was approximately constant when the damper was harmonically excited. The effect of temperature on the mechanical properties of VE dampers was analyzed by Chang et al. [10] where the empirical formulae describing the dependence of damper stiffness and the loss factor on temperature was established. Moreover, the self-heating process in VE dampers and layers was theoretically and experimentally investigated in [11][12][13][36][37][38][39][40].
The influence of temperature on the dynamic characteristics of structures with VE elements (dampers or layers) is very rarely explored. In [10], the dependence of the first natural frequency and non-dimensional damping ratio on different ambient temperatures was determined on the basis of experimental data. In [39], the effect of temperature on the dynamic characteristics of beams with a VE layer was investigated experimentally. The authors found the first few non-dimensional damping ratios to decrease with temperature increasing from 5 to 40 • C. Some information concerning the influence of temperature on the dynamic characteristics of structures with VE dampers was presented very recently in the papers [13,41]. Very few papers were also dedicated to the analysis of the influence of temperature on structures with dampers or layers in the time domain. In the Fig. 2 An illustration of the temperature-frequency principle paper [40], the results of the analysis of seismically excited structures in the time domain were presented and discussed. The finite element modeling temperature and frequency-dependent VE materials were considered in the paper [42]. In this paper, the frequency-temperature superposition principle introduced by Schwarzl and Staverman [43] is used to take into account the influence of temperature on the values of the parameters of rheological models adopted to the description of VE materials. A material to which this technique is applicable is said to be thermorheologically simple (see [43,44]). According to this principle, the curves of, for example, loss moduli versus frequency do not change their shape with temperature but are only shifted in the loss moduli-frequency plane. This is illustrated in Fig. 2.
It is assumed that the rod material is thermorheologically simple. According to the temperature-frequency principle, for s = iλ, the complex modules iλG(iλ), calculated for two different temperatures T and T 0 , where T 0 is the reference temperature, are equal. This can be written in the form of following equation: whereλ 0 is the reference frequency (see Fig. 2). The symbolα T is the shift factor. Equation (83) also results from the equality of ordinates in Fig. 2. According to relationship (83), the calculation of the function iλG(iλ) for the VE material working at a temperature T consists in shifting the argument of the function iλG(iλ) valid for the reference temperature T 0 in the following way (see Fig. 2 In a general case, two shift factors, i.e., the horizontal shift factor and the vertical shift factor, must be used for the description of some VE materials. As pointed out in the paper [36], for thermorheologically simple VE materials, only the horizontal shift factor is required. The horizontal shift factor is calculated from some empirical formulae. Very popular is the following William-Landel-Ferry (WLF) formula (compare [3]): where C 1 and C 2 are empirical constants and T = T − T 0 . Equation (85) is valid near the glass transition region (T g < T < T g + 100, according to the paper [45]). Moreover, it is assumed that above the temperature T g , the fractional free volume increases linearly with respect to time. Moreover, as the free volume of materials increases, its viscosity rapidly decreases (compare [46]). In the case of the fractional Zener model of the VE material, the temperature-frequency principle allows us to write: where Ē =Ē ∞ −Ē 0 and material parametersĒ 0 ,Ē ∞ ,τ andᾱ refer to the temperature T , whereas the parameters without the dashes refer to the reference temperature T 0 . Equation (86) is fulfilled when which leads to the conclusion that a change of temperature only causes a change of relaxation time in the following way:τ = τα T .
(88) Proceeding in a similar way, it is easy to find that relationship (88) is valid also for both the classic and the fractional Kelvin models. For the generalized Maxwell model, it is easy to obtain (for i = 1, 2, . . . , m) that Relationships (88) and (89) are crucial to the proposed approach because, in this simple way, it is possible to determine the parameters of VE materials for different temperatures.
It can be observed that temperature only changes the relaxation time (or times in the case of the generalized Maxwell model) and only the function θ(s) appearing in Eqs (22), (54) and (76) changes its values with temperature. However, this happens when the temperature changes are spatially homogenous for all structures. Moreover, the temperature influenced the frequency response function only through the function θ(s).

Results of representative calculations
The examples presented in this chapter have been selected so as to enable the authors to check the correctness of the developed software for frame structures built of VE materials. In addition, these considerations are subject to the influence of change in the parameters describing the VE material on the dynamic characteristics of systems.

Example 1: VE simply supported beam
The dynamic characteristics of a VE simply supported beam are presented in [47]. In this paper, the beam is considered as the non-local system but the results of calculation for the beam treated as the local system are also presented. In the case of the Kelvin-Voigt beam, the solution to Eq. (22) can be found analytically. Now, θ(s) = sτ and (22) has the form of the following quadratic equation: for which the solutions are These solutions can easily be compared with the analytical one, presented in [47] as Formula (29). In fact, both compared formulae are identical. The results presented in Table 1 agree very well for the simply supported beam and for the cantilever. However, the results presented for the fixed-fixed beam are different. Formula (91) is very simple. We check our calculation very carefully, and our final conclusion is that the results presented in [47] do not seem to be accurate for the fixed-fixed beam. It should be noted that in [47], the in-time solution is assumed in the form exp(iωt) while here we write exp(st) which justifies the differences in Table 1.  . 3 The ratios ω/ω e versus E ∞ /E 0 for the simply supported beam. ω 1 /ω e1 ( ), ω 2 /ω e2 ( ) Fig. 4 The non-dimensional damping ratios γ versus E ∞ /E 0 for the simply supported beam. γ 1 ( ), γ 2 ( )

Example 2: One-span VE beams
A family of representative calculation were obtained for a one-span beam of the length L = 4.0 m, the cross section dimensions b = h = 0.4 m and the mass per unit length m = 160.0 kg/m. The VE material of the beam is described by different models which will be specified later. The first results concern the beam made of the material modeled by the classic Zener model. The parameters of this model are: E 0 = 7.0 MPa, τ = 20.0 ms and the values of E ∞ change from E ∞ = E 0 up to E ∞ = 12E 0 . The first and second natural frequency ratios ω/ω e versus E ∞ /E 0 ratio are shown in Fig. 3. In Fig. 4, the dependence of the non-dimensional damping ratio of the first and second modes of vibration is shown versus the E ∞ /E 0 ratio. These results are for the simply supported beam. If that beam is treated as the elastic one with the elastic modulus E 0 = 7.0 MPa, then the first two natural frequencies are: ω e1 = 5.959 rad/s, ω e2 = 23.837 rad/s. The influence of the support conditions on is illustrated in Figs. 5, 6, 7 and 8. Figures 5 and 6 are for the fixed-fixed beam for which ω e1 = 13.509 rad/s, ω e2 = 37.238 rad/s while Figs. 7 and 8 are for the fixed-simply supported beam (ω e1 = 9.3096 rad/s, ω e2 = 30.169 rad/s. It can be seen that the type of support and the E ∞ /E 0 ratio have a significant influence on the dynamic characteristics of the beam. Basically, natural frequencies increase with the E ∞ /E 0 ratio but the progress of ω 2 /ω e2 ratio is faster than that of the ratio ω 1 /ω e1 . More complex is the change of non-dimensional damping ratios. Generally, γ 2 increase faster than γ 1 for E ∞ /E 0 < 3.5 but after that, the opposite tendency is observed. Moreover, γ 2 has its maximums for all the beam supports, whereas γ 1 constantly increase for the simply supported beam. The ratios ω/ω e versus E ∞ /E 0 for the fixed-fixed beam. ω 1 /ω e1 ( ), ω 2 /ω e2 ( ) Fig. 6 The non-dimensional damping ratios γ versus E ∞ /E 0 for the fixed-fixed beam. γ 1 ( ), γ 2 ( ) Fig. 7 The ratios ω/ω e versus E ∞ /E 0 for the fixed-simply supported beam ω 1 /ω e1 ( ), ω 2 /ω e2 ( ) The most complex results were obtained for the fixed-simply supported beam. In the range 9.385 ≤ E ∞ /E 0 ≤ 9.425, Eq. (22) has three real solutions instead of two complex conjugates and one real outside of this range. This means that the beam is overdamped, the motion of the beam is non-oscillating, and we cannot speak about Fig. 8 The non-dimensional damping ratios γ versus E ∞ /E 0 for the fixed-simply supported beam γ 1 ( ), γ 2 ( ) Fig. 9 The dependence of real solution of Eq. (22) on the ratio E ∞ /E 0 for the simply supported beam s 1 ( ), s 2 ( ), the fixed-fixed beam s 1 ( ), s 2 ( ) and the simply supported-fixed beam s 1 ( ), s 2 ( ) the natural frequency. Therefore, a jump of function ω 1 /ω e1 (E ∞ /E 0 ) can be seen in Fig. 8. All of this evidence leads to the conclusion that there are two bifurcation points in the s(E ∞ /E 0 ) curves drawn on the complex plane. The coordinates of the bifurcation points (denoted as s b and E ∞b ) could be determined analytically as well. This is presented above in Sect. 2. In the considered case of a fixed-simply supported beam, the following coordinates were obtained for bifurcation points which have the physical meaning: s b = −14.0972, E ∞ /E 0 = 9.3866 (point 1) and s b = −18.9945, E ∞ /E 0 = 9.4275 (point 2). These coordinates are in perfect agreement with the ones obtained numerically.
Equation (22) has also the real solution for the classic Zener model. In Fig. 9, the real solutions connected with two considered modes of vibration are shown in dependence on E ∞ /E 0 for all of the considered beam supports. For the simply supported-fixed beam, a jump of the s(E ∞ /E 0 ) curve is also observed. The values of s 1 for diversely supported beams, whereas the support conditions have a rather minor influence on the s 2 values.
The influence of the fractional derivative α on the natural frequencies and the non-dimensional damping ratios is also examined. The results of the calculations are presented in Figs. 10, 11, 12 and 13. In Fig. 10, the first and second natural frequency ratios ω/ω e are shown in dependence on the order of the fractional derivative α. The results for the simply supported beam are drawn as the solid lines, for the fixed-simply supported beam as the solid line with small crosses, and those for the fixed-fixed beam are presented as the solid line with small circles. In a similar way, Fig. 11 shows the two first non-dimensional damping ratios. The above-mentioned results are for E ∞ /E 0 = 10/7, and those for E ∞ /E 0 = 5 are in Figs. 12 and 13, presented in a similar way. For E ∞ /E 0 = 10/7, the natural frequencies decrease with increasing values of α for all considered types of Fig. 10 The ratio ω/ω e versus the order of fractional derivative α for different beams and the ratio E ∞ /E 0 = 10/7. The simply supported beam ω 1 /ω e1 ( ), ω 2 /ω e2 ( ), the fixed-simply supported beam ω 1 /ω e1 ( ), ω 2 /ω e2 ( ) and the fixed-fixed beam ω 1 /ω e1 ( ), ω 2 /ω e2 ( ) Fig. 11 The non-dimensional damping ratios γ versus the order of the fractional derivative α for different beams and the ratio E ∞ /E 0 = 10/7. The simply supported beam γ 1 ( ), γ 2 ( ), the fixed-simply supported beam γ 1 ( ), γ 2 ( ) and the fixed-fixed beam γ 1 ( ), γ 2 ( ) beam supports but changes are rather minor. In contrast to this, for E ∞ /E 0 = 5, the first natural frequency also decreases with increasing α but the second natural frequency increase with α. Moreover, these changes are much more pronounced. That natural frequencies increase with increasing α is justified by the fact that the elastic properties of the VE material decrease with increasing values of α.

Fig. 16
The ratios ω/ω e versus T for the fixed-simply supported beam ω 1 /ω e1 ( ), ω 2 /ω e2 ( ). Moreover, E ∞ /E o = 3,α = 0.7, τ = 20 ms, C 1 = 9.23, C 2 = 141.2 K, T 0 = 20.0 • C 7.1 Example 3: The fixed-simply supported beam-influence of temperature The influence of temperature on the first and second natural frequencies and the non-dimensional damping ratio is presented in Figs. 16 and 17, respectively. The results are shown for the fixed-simply supported beam for which: E ∞ /E o = 3.0, α = 0.7 and τ = 0.20 ms is chosen as data at the reference temperature T o = 20 • C. Moreover, C 1 = 9.23 and C 2 = 141.2 K. The solid lines present the results for the first natural frequency and the first non-dimensional damping ratio, whereas the dashed lines show the results for the second ones.
Both the natural frequencies and non-dimensional damping ratios can change their values over a wide range. The natural frequencies significantly decrease with increasing temperatures. The most important changes are observed for temperature values between − 20 and 35 • C. Outside this range, the natural frequencies are nearly constant. For higher temperatures, the natural frequencies are nearly equal to the natural frequencies of the corresponding elastic beam and approximately 75% higher than those frequencies of the elastic beam for temperatures lower than −10 • C. The non-dimensional damping ratios are also strongly influenced by temperature but these functions have their maximums at about 15 • C and 21 • C for the first and second damping ratios, respectively. Outside the temperature range (−15 • C, 50 • C), the damping ratios are very low and not greater than 0.005. On the other hand, the maximums of damping ratios are very high at approximately 0.37. The presented results illustrate how the dynamic characteristics of the beam can change with temperature.

Concluding remarks
The new method for determining the dynamic characteristics of frame structures made of materials exhibiting VE properties is presented. The method can be used when the VE material is described by different rheological models (both classic and fractional ones). However, the structures must be homogenous and one rheological model must be able to describe the material of the frame. This means that the proposed method is not general but can be applied to an analysis of a very important class of structures. The dynamic characteristics of the rod treated as the continuous system are determined. Moreover, the ordinary and dynamic finite element formulations for the VE structures of planar frame are presented. The influence of temperature on the dynamic properties of structures can be easily taken into account using the temperature-frequency principle.
Some interesting dynamic properties of the considered VE structures have been discovered. It was found that the modes of vibrations are real ones. Moreover, it was found that, for some rheological models of the frame material, both the viscous and non-viscous though always real modes of vibration exist. These modes are identical if they come from one solution to a linear eigenvalue problem defined for the corresponding elastic structure. Moreover, a specific form of orthogonality conditions was found which is valid for the considered type of VE structures. It is also valid for frames made of a VE material described with the help of fractional derivatives. The above-mentioned orthogonality conditions enable the calculation of the frequency response functions in a very efficient way.
The method is very attractive from the computational point of view. In contrast to the general case when a highly time-consuming solution to the nonlinear eigenvalue problem is needed (see, for example, [13,15,16,18,20,25,26,29]), the proposed method requires a solution to the linear eigenvalue problem for the corresponding elastic structure and a solution to one nonlinear algebraic equation for each needed nonlinear eigenvalue. In contrast to the presented method, the solution to a large system of nonlinear equations is needed when typical methods, such as the continuation method or the perturbation method, are used.
The presented results of calculation illustrate how the dynamic characteristics of beams depend on many factors: the support conditions, modulus ratio, order of fractional derivative and relaxation time.
For all of the considered boundary conditions, the non-dimensional damping ratios have their maxima for the specific E ∞ /E o ratio. Except the first frequency ratio for the simply supported beam, the other frequency ratios grow with the increasing E ∞ /E o ratio.
The frequency ratio ω/ω e depends on the order of fractional derivative α in a different way. For the relatively low E ∞ /E o ratio, natural frequencies slightly decrease when α increases. In contrast to this, for high E ∞ /E o ratios, the first natural frequency significantly decreases while the second natural frequency considerably increases with increasing α. The non-dimensional damping ratio increases with α in most cases but the changes are more significant for higher E ∞ /E o ratio.
The frequency ratios grow with relaxation times. The corresponding non-dimensional damping ratios also increase for low E ∞ /E o ratio and have their maxima in the case of high E ∞ /E o ratios. Acknowledgements This study has been supported by the Poznan University of Technology, Institute of Structural Engineering, carried out in the year 2019. This support is gratefully acknowledged.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http:// creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.

Appendix A
A brief description of the most popular rheological models used as the models of the VE material is presented below, where the constitutive relationships in the time and the Laplace domain were written. For the classic Kelvin model, the above-mentioned constitutive relationship in the time domain is where τ = η/E 0 and η is the relaxation time and the viscosity, respectively. In the Laplace domain which means that The derivative of function θ(s) = sτ is θ(s) = sτ and, in that case, the function f (s) defined in (22) is a holomorphic function.
For the fractional Zener model, the analogous relationships are: from which it is easy to obtain the followinḡ For this model It is easy to verify that, for all of the considered models, is it possible to write sḠ(s) = (1 + θ(s))E 0 . (A.20)

Appendix B
In the local coordinate system, the elemental stiffness matrix of the rod executing both axial and transverse vibration is in the following form: The numerical procedure used to solve the nonlinear eigenvalue problem requires the calculation of the dynamic stiffness matrix K of the elastic structure. The structure can be modeled using finite elements with different masses per unit length m e and the geometric parameters l e and J e . In consequence, λ e could be different for each finite element. Therefore, the reference finite element is introduced with the parametersm r ,l r ,J r and with the reference characteristic value λ related to these parameters in the following way: