Fixed-Point Iteration Method for Uncertain Parameters in Dynamic Response of Systems with Viscoelastic Elements

: The paper presents a method for determining the dynamic response of systems containing viscoelastic damping elements with uncertain design parameters. A viscoelastic material is characterized using classical and fractional rheological models. The assumption is made that the lower and upper bounds of the uncertain parameters are known and represented as interval values, which are then subjected to interval arithmetic operations. The equations of motion are transformed into the frequency domain using Laplace transformation. To evaluate the uncertain dynamic response, the frequency response function is determined by transforming the equations of motion into a system of linear interval equations. Nevertheless, direct interval arithmetic often leads to significant overestimation. To address this issue, this paper employs the element-by-element technique along with a specific transformation to minimize redundancy. The system of interval equations obtained is solved iteratively using the fixed-point iteration method. As demonstrated in the examples, this method, which combines the iterative solving of interval equations with the proposed technique of equation formulation, enables a solution to be found rapidly and significantly reduces overestimation. Notably, this approach has been applied to systems containing viscoelastic elements for the ﬁ rst time. Additionally, the proposed notation accommodates both parallel and series con ﬁ gurations of damping elements and springs within rheological models.


Introduction
The process of designing and analyzing structures involves the consideration of a wide range of factors that can influence their behavior under dynamic conditions.Uncertainties concerning design parameters, materials, production processes, and environmental conditions are significant factors that can greatly impact the safety, reliability, and performance of structures.Traditional approaches to dynamic structural analysis often overlook these uncertainties associated with design parameters, leading to inaccuracies in predicting structural behavior under real operating conditions.Therefore, addressing uncertainties is a crucial aspect of modern design processes.
There are some main approaches to incorporating uncertainties associated with design parameters.The first is the probabilistic approach, which assumes that uncertainties in design parameters are modeled using statistical parameters.This allows for the determination of the probability of various scenarios, and examples of its application can be found in works by Ghanem and Spanos [1], Capillon et al. [2], as well as Wang et al. [3], including research within the context of structures with viscoelastic elements.However, many cases present challenges such as insufficient statistical data or difficulties in precisely determining probability distributions for all design parameters.In such cases, the second approach, involving interval analysis, is utilized.Here, design parameters are defined by interval values.This article employs interval analysis to investigate the impact of uncertainties in design parameters on the dynamic response of structures.
One of the simplest tools that can be applied is the Monte Carlo method and the vertex method.The Monte Carlo method involves repeatedly sampling from specified intervals for individual design parameters and then conducting an analysis for each sample (Jiang et al. [4]).While effective, its main drawback is the high computational cost, especially for complex construction models and a large number of parameters.On the other hand, the vertex method relies on representing intervals using vertices and finding solutions for all their combinations (Qiu et al. in [5,6]).This method is also time-consuming, which is why, in recent decades, efforts have been focused on finding other methods to reduce the computational burden.
Chen et al. in [7,8] and Yang in [9] presented a method for computing uncertain eigenvalues of systems with interval parameters.They proposed the utilization of matrix perturbation theory and interval extension theory.Both undamped and damped cases were considered.A similar approach was applied in [10] by Qiu and Wang but for determining the dynamic response of structures.The authors also conducted a comparison of the obtained results with the probabilistic approach, based on both mathematical proofs and numerical simulations.The range of obtained results encompassed those obtained by the probabilistic method.
An alternative to perturbational methods has been presented in [11] by Xia et al.The authors introduced a method based on the vertex method, which was applied to the first-order deviation of dynamic response.They demonstrated that this proposed method reduces overestimation and is also more accurate than perturbational methods.However, it is worth noting that second-order components in the equations of motion were omitted in their considerations; hence, the method may not be suitable when parameter variability is high.
Sim et al. in [12] and Moens and Vandepitte in [13] applied modal interval analysis, combining modal analysis and interval calculus.The authors presented a method for finding eigenfrequencies, modal shapes, and frequency response functions (FRFs).
Yaowen et al. in [14] determined the frequency response function for systems with uncertain parameters using the Laplace transform.This transformation converted dynamic equilibrium equations into a linear system of algebraic equations.The authors also applied the element-by-element technique, which significantly reduced overestimation resulting from the dependencies between interval parameters.However, this method is only suitable for handling small uncertainties.
The procedure for deriving approximate explicit values of FRFs for systems with uncertainties was presented in [15] by Muscolino et al.The effectiveness of the Rational Series Expansion in estimating the lower and upper bounds of FRF values was demonstrated.The proposed method also effectively mitigates overestimations.
In the work [16] written by Wei et al., an approach based on two-dimensional Chebyshev polynomials was introduced for estimating the dynamic response of structures with interval uncertainties.The authors utilized the decomposition of the dynamic response function into univariate and bivariate functions and fitted them using Chebyshev polynomials.This approach facilitates the transformation of the dynamic analysis of nonlinear systems with uncertain parameters into solving simpler one-dimensional and two-dimensional equations.
In studies concerning uncertain design parameters, obtaining the response of structures, particularly in dynamic problems, is greatly facilitated by expanding the sought value into a Taylor series around the central value.This approach utilizes the sensitivities of the sought function with respect to changes in uncertain design parameters.In the work [17], Qiu et al. applied this approach to estimate the dynamic response of damped systems.In reference [18] written by Chen et al., it was shown that using a second-order Taylor series enables the accurate estimation of structural responses even with large uncertainties, although the focus was on undamped systems.In reference [19], Fujita and Takewaki proposed the identification of critical combinations of uncertain parameters by maximizing the objective function through second-order Taylor series expansion.The structure under consideration included dampers, with uncertain parameters extending to damper characteristics as well.Additionally, in the works [20] by Łasecka-Plura and Lewandowski and [21] by Łasecka-Plura, structures featuring viscoelastic elements were analyzed.These studies presented algorithms that notably reduced the number of combinations needed to determine the lower and upper bounds of the sought dynamic response of the structure.
One method that enables the incorporation of uncertainties in design parameters is stochastic analysis (Kamiński et al. in [22,23]).This approach, particularly when parameters are described as interval numbers, was applied by Muscolino and Sofi in [24] using affine arithmetic.This application helped reduce the overestimation of the determined structural response.Instead, Wang et al.,in [25], introduced a hybrid approach where the boundaries of probabilistic characteristics are determined through interval arithmetic.Additionally, an approach utilizing Chebyshev polynomials to solve problems with multiple uncertain variables was introduced.
The authors also address scenarios where the analyzed systems are described using mixed random and interval parameters.In reference [26] by Zhao, a method was proposed to determine the lower and upper bounds of the mean value and standard deviation assuming the presence of mixed uncertain parameters.This method is based on expanding a series of random and interval quantities with respect to uncertain parameters.A comparison with the Monte Carlo method was also conducted, demonstrating a significant reduction in computational costs.In reference [27], Feng et al. considered the time response of structures and proposed a hybrid model combining generalized polynomial chaos theory and the Legendre metamodel method.On the other hand, Zhao et al., in refence [28], presented a method based on genetic algorithms for computing optimization problems with hybrid uncertain parameters.
Interval analysis finds broad application when dealing with uncertain design parameters, extending to other problems related to the analysis of structural behavior under dynamic loads.In [29] by Sofi et al., it was applied to analyze the reliability of systems with viscoelastic dampers; in [30] by Wang et al., it was applied to assess the reliability of systems with active vibration control; and in [31], Fang and Huang utilized it in the analysis of structural damage detection involving various types of uncertainties.
Uncertainty consideration is also discussed in the context of optimization.In [32], Wang et al. proposed a method for manipulator trajectory control, taking into account uncertainties.This method utilizes a non-probabilistic uncertainty quantification approach, ensuring precision and effectiveness in evaluating system response.On the other hand, in [33], Wang et al. presented topology optimization for continuous structures, considering uncertainties and time-dependent reliability.It employs a parameterized level set method based on reliability criteria.Meanwhile, in [34], Liu and Wang introduced a hybrid optimization method for vibration controller design, considering mixed uncertainties.Uncertainties are treated as intervals and fuzzy models.
As part of the analysis of systems with uncertain parameters, the issue of interval sensitivity is also considered, where the influence of individual uncertain input parameters on the range of obtained results is determined.In [35], Moens and Vandepitte demonstrated a procedure for determining the sensitivity of frequency response function bounds.
This article analyzes viscoelastically damped systems with uncertain design parameters.Methods for incorporating uncertainties in such systems, as described in the literature, can be divided into probabilistic methods, methods based on fuzzy theory, and methods where uncertainties are modeled as interval numbers.The first group includes works such as [36,37] written by Hernández et al., which analyze dynamic characteristics for systems with uncertainties, and [3,38] by Wang et al. and [2] by Capillon et al., where the authors consider random design parameters in determining frequency response function.Fuzzy theory, on the other hand, is applied in works [39,40] written by Nasab et al., where frames with viscoelastic dampers were analyzed.The third approach is described, among others, by Fujita and Takewaki in [19], Łasecka-Plura and Lewandowski in [20], Łasecka-Plura in [21], and Sofi et al. in [29], which will be discussed later in this section.A detailed description of incorporating uncertainties in structures with viscoelastic elements can be found in the review paper [41] by Lewandowski et al.The choice of the appropriate approach depends on the description of uncertainty.In the case of viscoelastic materials, it is often difficult to determine the parameter distribution, making interval analysis a good solution as it requires knowledge of only the lower and upper bounds of parameter variability.
In this article, the dynamic response of structures with viscoelastic elements is analyzed.It is assumed that the behavior of viscoelastic material is described using rheological models, both classical and fractional.The method presented considers uncertainties in the design parameters of the structure as well as damper parameters.A method for determining the frequency response function is presented, which is one of the most important tools for assessing the dynamic response of structures.However, there are very few publications addressing this issue for structures with viscoelastic elements.Interval analysis is employed in this article, addressing a drawback of this approach, namely significant overestimation.To avoid this, the element-by-element technique is applied.Previous works, ref. [42] by Modares et al., ref. [43] by Muhanna and Zhang,and [44] by Yaowen et al., have shown that this technique reduces overestimation in both statics and dynamics problems.Its application to systems with viscoelastic elements is presented here for the first time.The formulation of equations of motion using the element-by-element technique is general and applicable to various rheological models.Laplace transformation is used to express the equations of motion in the frequency domain, representing them as a system of interval equations.Their solution involves individual elements of the FRF matrix, solved using Brower's fixed-point theorem.Three methods are compared in computations: the proposed method, the vertex method, and the direct arithmetic of interval numbers.Direct analysis leads to significant overestimations, while the proposed method aligns more closely with the vertex method and reduces computational costs significantly.The examples presented demonstrate that the proposed method is a valuable tool for assessing the lower and upper bounds of the response of structures with viscoelastic dampers.Additionally, the obtained results in one example were compared with those obtained using Taylor series expansion.
A novelty of this work is the application of the fixed-point iteration method to find the lower and upper bounds of the frequency response function for systems with viscoelastic damping elements, also described using fractional derivatives.The presented method can be applied to account for the influence of uncertainties in the parameters of viscoelastic damping devices on the response of the device or the entire structure, which is an advantage compared to traditional methods that do not consider parameter uncertainties.
The article is organized as follows: Section 2 provides a brief introduction to interval analysis; Section 3 presents the formulation of equations of motion for systems with viscoelastic elements, considering uncertain design parameters; Section 4 demonstrates the application of Brower's fixed-point theorem to solve systems of interval equations for determining the frequency response function; Section 5 presents the method of overestimation reduction, involving the formulation of equations of motion using the element-by-element technique with the use of Lagrange multipliers; Section 6 presents examples; and Section 7 concludes the article.

Introduction to Interval Analysis
In this article, uncertainties in design parameters are represented as interval numbers.A comprehensive discussion of interval arithmetic basics is available in [45] by Neumaier and in [46] by Moore.This section will focus on presenting the most important fundamental information.
The design parameters of the considered system can form a set, which can be denoted as a column vector  col  ,  , … ,  , where s is the number of parameters.The values of these parameters are not precisely determined, but the lower and upper bounds of their variability are known, and they can be written as interval numbers.An uncertain parameter  can be represented as follows: where • denotes an interval number,  represents the lower bound of the design parameter, and  represents the upper bound of the design parameter.It can also be expressed in the following way: or where    2 ⁄ denotes the central value of the design parameter,    2 ⁄ represents half of the interval determining the parameter, and  stands for the uncertainty of the parameter  .In this article, the notation specified by Equation ( 3) is employed.Basic operations on interval numbers are presented in Appendix A.
In general, computations involving interval numbers are significantly more complex and time-consuming compared to real numbers.This complexity arises from the need to consider all possible combinations of lower and upper bounds of parameters, followed by selecting the lower and upper bounds of the resulting value.Consequently, this leads to increased computational overhead.
An additional challenge in analyzing engineering problems using interval analysis is the potential for significant overestimation resulting from the direct application of interval arithmetic.This occurs because, in formulating equilibrium equations for complex systems, one physical quantity may sometimes be represented by more than one interval number.Additionally, uncertainties of interval numbers may overlap, resulting in a much broader interval in the computation outcome than the one containing all possible solutions.Therefore, reducing overestimation is a key issue to address when employing interval arithmetic.This problem is discussed, among others, by Neumaier in [45] and by Rump in [47].In this article, overestimation reduction was achieved through the application of the element-byelement technique and by transforming the solution into an iterative form, which will be discussed in detail based on the damping parameters of damping elements in Section 5.

Dynamic Response of Systems with Viscoelastic Elements
The analysis will focus on a structure with viscoelastic damping elements, described using rheological models that incorporate fractional derivatives.The general equation of motion can be expressed as follows: where  and  are the mass and stiffness matrices, respectively,  is the damping matrix associated with the structure, and    is a matrix of damping kernel functions whose form depends on the adopted model.  represents the vector of excitation forces and   is the vector of displacements.In structural dynamics for matrix    , rheological models are commonly used, with fractional models being employed in this article.As part of this analysis, the Scott-Blair viscoelastic element (Figure 1) will be introduced.Its behavior is described by the following equation: where  and  are the parameters of the described element 0  1 , the symbol  • denotes the fractional derivative of order  with respect to time, and   is the difference in displacements of the damper nodes.  denotes the force in the viscoelastic element.A more extensive discussion of fractional calculus can be found in [48].The following Riemann-Liouville definition of the fractional derivative was applied in the article: where Γ denotes the gamma special function.As demonstrated in [49][50][51], the application of this definition is justified in the context of viscoelastic materials.In this article, Kelvin or Maxwell models will be considered, both classical and fractional (Figure 2).Classical models are described in [52] by Lewandowski et   The force in the damping element, in the case of the Kelvin model, is described by the following equation: and for the Maxwell model, it is described as follows: where  and  are the damping coefficients of the Kelvin and Maxwell model, respectively, and  and  are the stiffness parameters of the Kelvin and Maxwell model, respectively.Additionally,    ⁄ , and by achieving α = 1, one can obtain the classical models.
After applying the Laplace transformation, Equation ( 4) can be written as follows: where   and   represent the Laplace transform of   and   , respectively, and   is the matrix of damping kernel functions in the Laplace domain.Its form depends on the adopted rheological model of the viscoelastic element, the considered structure, and the placement of the viscoelastic element.It can be expressed as the sum of matrices for individual viscoelastic elements, denoted as   ∑   (where r denotes the number of viscoelastic elements).
Considering steady-state harmonic responses, the harmonic external forces can be described as follows: where  is the vector of excitation force amplitudes, and  is the frequency of excitation.
The displacement of the structure can then be expressed as follows: Substituting ( 10) and (11) into Equation ( 4), we obtain the following formula: where   is the dynamic stiffness matrix, which can also be obtained by substituting into Equation (9)  i.

𝐃 𝜆 𝜆 𝐌 i𝜆𝐂 𝐊 𝐆 𝜆
We assume that all design parameters can vary within certain specified bounds; therefore, we can represent them as interval numbers.The equation of motion in general form can be written as follows: or where The solution to Equation ( 15) can be expressed as follows: where   denotes the interval frequency response function and can be determined as follows: The individual elements of the matrix   with indices i and j can be interpreted as the frequency response function determined for the degree of freedom with index i when unit harmonic excitation acts in the direction of the dynamic degree of freedom of the structure with index j.Determining the FRF matrix from Equation ( 18) requires the direct application of interval arithmetic.This method will later be referred to as the direct interval method (DIM).However, as demonstrated in the examples presented, this approach leads to significant overestimations, as described by Yaowen et al. in [44] and Moens and Vandepitte in [13].
To mitigate the overestimation, a method will be applied where the calculation of the FRF matrix will be divided into steps, computing each column of the matrix   sequentially.Since the dynamic response of the structure, understood as vector   , can be obtained using Equation (17), it is worth noting that if we represent the vector   as the first column of the matrix   will be determined.Subsequent columns can be obtained by substituting value 1 at successive positions in vector .This approach involves solving interval equation systems to compute the FRF matrix.Various iterative methods that can be applied for this purpose are presented by Rump in [47,56] and Jansson in [57].
The specific method used in this work will be described in the next section.

Interval Equations System Solving
In the paper, Brower's fixed-point theorem was applied to solve the system of linear interval equations (Rump [47], Neumaier [58]): after substituting the appropriate vector   .Equation ( 20) can be transformed into a fixed-point equation in the following form: where and  is a nonsingular matrix.In Equation ( 22), all quantities except  are interval numbers.Using Brower's fixed-point theorem, we can express the following: The solution to Equation ( 23) also satisfies the original Equation (20).Therefore, it is possible to write the following (Rump [59]): where   denotes interior of  .
If we denote that    , where  is the computed vector for the central values of the parameters, and  is the uncertainty of the solution, then Equation ( 24) can be expressed as follows (Neumaier [45]): The equation used for the iterative determination of the uncertainty of the solution can thus be written as follows: According to work [45] by Neumaier, the matrix  is defined as    .The stopping criterion for the iteration is as follows: The final value of the solution is obtained based on the following equation: The above algorithm can be successfully applied to solve systems of linear interval equations, as demonstrated in [43] by Muhanna and Zhang.

Application of the Presented Method to Multi-Degrees-of-Freedom Systems with Viscoelastic Elements
The matrix-form equation of motion for systems with viscoelastic elements is given in the form (12).However, such representation leads to interval quantities associated with one physical quantity appearing more than once, which can cause significant overestimations.As demonstrated Muhanna and Zhang in [43] for statics and Yaowen et al. in [44] for dynamics, overestimation is reduced by applying the element-by-element technique.This approach is detailed in [44] by Yaowen et al.As a result of employing this technique, we obtain a singular stiffness matrix.Hence, it is necessary to introduce an additional constraint matrix  for the considered system and apply the Lagrange multiplier method  .Then, the equation of motion looks as follows: where interval vectors and matrices are constructed using the element-by-element technique.This technique will be demonstrated based on the construction of the component matrix   .Let us consider a three-story frame, as shown in Figure 3, with two dampers described by the Kelvin model.The schematic illustrates the forces  , representing the interactions of individual dampers with each floor.The matrix containing uncertain damping parameters  will be considered.According to [51] by Lewandowski and Pawlak and [52] by Lewandowski et al., for one damper k, the matrix can be expressed as follows: where  is the uncertain damping parameter of damper k,  is the central value of parameter  , and  denotes the uncertainty of parameter  .
If we consider two dampers, it is possible to express the matrix  in the following form: It is worth noting that Formula (31) aims to illustrate only the idea of the element-by-element technique using the example of coefficients for two dampers, and it is not the final form of the matrix appearing in the equation of motion.From the above expression, it is clear that one uncertain parameter associated with the damper appears four times, meaning the uncertainties  overlap and cause overestimation.According to the above notation, the remaining matrices can be expressed as follows: To reduce overestimation when multiplying any matrix by a displacement vector, the following transformation will be applied, expressed in the general form for quantities  and  : The above transformation ensures that each uncertain parameter associated with quantity i appears only once.Equation ( 33) can be rewritten in the following form: or in a more compact form, as follows: where  ∑   ,   ,  i ,    and    .Assuming that Equation (36) can be treated as a system of linear interval equations, Equation ( 26) can be written in the following form: and after considering (36) and   , we obtain the iterative equation in the following form: In Equation (38), in the case of multiplication   , we apply transformation (34).

Single-DOF System: Kelvin Model of Viscoelastic Element
To illustrate the method's idea, we will first discuss the single-degree-of-freedom system with a Kelvin damper, both fractional and classical, as shown in Figure 4.The equation of motion (15), assuming that the system parameters are described using interval numbers, now takes the following form: where  and  are understood as the mass and stiffness of the system, respectively, and  ,  , and  represent the stiffness, damping coefficient, and the order of the fractional derivative describing the parameters of the Kelvin damper, respectively.  denotes the forcing function, and   represents the displacement of the considered system.The task is to find the system response in the form of displacement values.We assumed the following parameters for the analyzed system:  40 kg ,  1000 N/m,  500 N/m.Various cases of damping coefficients were analyzed, with the data shown in Table 1.It was assumed that the forcing magnitude  1 N and is not affected by uncertainties.Only the stiffness and damping parameters of the viscoelastic element are uncertain, with 5% uncertainty.Table 2 presents the parameter values in the form of interval numbers.The response of the structure was calculated using the vertex method (red line), DIM (green line), and the method described in the article (blue line).The results are presented in Figures 5-8.The zoomed-in figures show only the region near the resonance zone because in other areas, the results are practically overlapping.
In the first two cases, the classical Kelvin model was considered (Figures 5 and 6), while in the next two, the fractional model was analyzed (Figures 7 and 8).Im(q()) Re(q()) direct interval method vertex method presented method The results for the presented method and for the DIM almost overlap, especially when the system's damping is higher.The obtained results are also very close to those of the vertex method, indicating that the presented method can assess the influence of uncertain design parameters on the system's response.all cases, the results of the presented method lie outside the interval obtained from the vertex method, indicating that they are on the safe side.The example shown demonstrates that the presented method can be applied to determine the response of a system with a viscoelastic element described using both classical and fractional Kelvin models.

Single-DOF System-Maxwell Model of Viscoelastic Element
The system under consideration is a single-degree-of-freedom system, as depicted in Figure 9, where the viscoelastic element is described using the Maxwell model.The presented method requires the equations of motion to be formulated in a way that describes the relationship between the central value of the parameter and its uncertainties, as outlined in Equation (31).In the case of the Kelvin model, as discussed in the previous section, achieving this is relatively straightforward.However, for the Maxwell model, it is necessary to use a representation involving an internal variable indicated in Figure 9 as  (Lewandowski et al. [52]).The equations of motion then take on the following form: or the following matrix form: From this expression, it is evident that even in a simple case like a single-degree-of-freedom system, when dealing with the damping parameter of the Maxwell model, there might be overestimation.This results from the fact that the interval number describing one parameter appears four times in the expression (refer to the second term of Equation ( 41)).To investigate the potential for reducing overestimation with the applied method, two cases were examined in the computed examples.In the first case, the transformation described by Equations (34) was not applied (Figures 10 and 11, while in the second case, the transformations (34) were applied (Figures 12 and 13).The analyzed system is characterized by the following parameters:  40 kg,  1000 N/m, and  500 N/m.Two different damping parameter values were assumed, as detailed in Table 3.Similar to example 6.1, the parameters of the viscoelastic damper have 5% uncertainty.The parameter values are listed in Table 4.     Im(q()) Re(q()) direct interval method vertex method presented method The results are calculated using three methods: vertex, DIM, and the presented method.
In the cases considered without ransformation (34), the results from the presented method and the DIM are almost identical, leading to overlapping plots in most instances.However, it is important to note that there is an overestimation, particularly evident in the example shown in Figure 10.In this case, the plots do not resemble the results from the vertex method at all, indicating complete inaccuracy.On the other hand, applying transformation (34) brings the results obtained using the fixed-point iteration method closer to those from the vertex method, leading to a significant reduction in overestimation.
Both in examples 6.1 and 6.2, the results are magnified for the resonance zone, as the outcomes in this zone are significant from an engineering perspective.The differences visible in the plots are thus emphasized through magnification.Based on the presented examples, it is also noticeable that the presented method provides more accurate results when the nondimensional damping ratio of the analyzed system is higher.

Frame with Kelvin Models
The structure under consideration is a frame with dampers described by the fractional Kelvin model, as shown in Figure 14.The parameters of the frame are as follows: the mass of the floor is  1000 kg , and the stiffness of the storey is  10000 N/m .Structural damping is assumed to be proportional and described by the relationship  0.001 1.0.Additionally, all dampers are assumed to be the same and described by the following parameters:  1000 N/m ,  500 Ns /m , and  0.8 .The natural frequencies and the nondimensional damping ratios for the analyzed frame are presented in Table 5.The parameters of the dampers  and  are uncertain, with 5% uncertainty assumed for each parameter independently.The parameters of the dampers can be described as follows: where  0.05, 0.05 .The matrix  appearing in Equation ( 36) has the following form: where the matrices  and  contain the central values of the damping coefficients and stiffness coefficients of the dampers, respectively, while  and  contain the uncertainties of these parameters.All matrices in Equation ( 43) are constructed using the element-byelement technique, and transformation (34) has been applied in the calculations.The real and imaginary parts of selected elements of the FRF matrix are shown in Figures 15-18.Based on Figures 15-18, it is evident that the presented method yields results close to the vertex method.Moreover, compared to the DIM, it significantly reduces overestimation resulting from the direct analysis in interval arithmetic, especially at lower frequencies.This confirms the advantages of the presented method.Additionally, considering the computational cost of the presented method and the vertex method, the presented method also holds an advantage.In the case of the vertex method and six varying parameters, the number of cases to be considered for each analyzed frequency is 2 64, whereas with the presented method, the solution is typically obtained after a few iterations and is relatively straightforward to implement.In Tables 6 and 7, the relative error for the presented method and the direct method is shown for selected values of frequency and two uncertainty values: 5% and 1%.The absolute errors of the magnitude for  were calculated.The vertex method was adopted as a comparative method.Im(H44()) Based on the results presented in Tables 6 and 7, it can be observed that the error values are always smaller in the case of the presented method.The largest differences are for lower frequencies, which confirms the observations from the  plots, where the largest overestimations for the direct method are indeed in this area.The obtained results confirm the superiority of the presented approach over the direct interval method.
In Table 8, the computation time for H11 is presented for various uncertainty values.Based on the example provided, it can be observed that the discussed iterative method is more efficient, the smaller the uncertainties are.This is associated with a lower number of iterations needed to find a solution.On the other hand, the vertex method depends on the number of parameters, and its computational cost is similar for each uncertainty case.It can also be noticed that the presented method is more advantageous in terms of computation time.The calculations were performed using a laptop with an Intel Core i7-7500U CPU 2.70 GHz and 16 GB RAM.The main aim of this article is to present a method allowing for the consideration of uncertainties in the parameters of viscoelastic damping elements.However, the presented method is general and can also account for uncertainties in structural parameters.There are no limitations regarding the number of uncertain parameters; however, the convergence and accuracy of the method decrease as the number of parameters increases.At the same time, the efficiency of the comparative method, such as the vertex method, also significantly decreases as the number of uncertain parameters increases, so despite some inconveniences, the presented approach remains advantageous.To demonstrate that it is possible to account for different numbers of parameters and uncertainties, the same frame is considered, but in addition to the uncertainties in dampers parameters, the stiffness of the stories and the mass of the floors are also uncertain.The stiffness of each story and the mass of each floor can change independently, so together with the uncertain damper parameters, the number of changing parameters is 14.It is worth noting that in the case of the vertex method, the number of combinations for each considered frequency is now 2 16,384.It was assumed that the uncertainty of the damper parameter is 1%, the uncertainty of the stiffness of the stories is 0.5%, and the uncertainty of the mass of the floors is 0.75%.The intervals of the design parameter values are presented in Table 9.The obtained results of  are shown in Figure 19.Therefore, it can be observed that the method can also be applied for the prediction of the dynamic response of structures when the number of parameters is large or when the assumed variabilities are different.The results obtained using the presented method were also compared with those obtained by a method based on Taylor series expansion.Predicted values of the frequency response function were computed for the case of considering the first two of the Taylor series (Łasecka-Plura and Lewandowski [20]).This approach is briefly discussed in Appendix B. It involves calculating the sensitivity of the frequency response function with respect to changes in successive parameters and then selecting the largest and smallest predicted increment.In Figures 20-23, a comparison between the vertex method, the presented method, and the Taylor series expansion-based method (black line) is depicted.Comparative calculations were conducted for two damping parameter variabilities of 1% (Figure 20) and 2% (Figure 23).Taylor-based method vertex method presented method Im(H11())  In the first case, the results obtained from the presented method and the Taylor series expansion-based method are very similar, indicating that both provide similar approximations.Figures 21 and 22 show enlarged parts of the graph.It can be observed that the presented method yields results closer to the vertex method than the Taylor series expansion-based method.In the case of a parameter variability of 2%, the Taylor series expansion approach for higher frequencies yields results deviating from the vertex method, unlike the presented method.Therefore, the presented method has an advantage over the Taylor series expansion-based method when considering the first two terms of the series.

Frame with Maxwell Model
The frame shown in Figure 14 is considered, but with dampers described using the Maxwell model.The parameters of the frame remain the same as in example 6.3, and the structural damping is assumed in a similar manner.The dampers are characterized by the following parameters:  1000 N/m ,  500 Ns /m , and  0.8 .Table 10 presents the natural frequencies and nondimensional damping ratios for the analyzed frame.As previously mentioned, the focus is on the influence of uncertainties in the parameters of viscoelastic elements.Therefore, it is assumed that the parameters and the dampers  and  are uncertain, with uncertainty of 5%.It is further assumed that these parameters can change independently and are described as follows: where  0.05, 0.05 .This time, unlike the Kelvin model, to express the relationship between uncertainty and the central value of the parameter in Equation ( 36) as a multiplication, it is necessary to use the formulation of equations of motion with internal variables.A comprehensive explanation of constructing such equations of motion for a frame with viscoelastic dampers can be found in [52]       In this case, it can also be observed that the applied method closely resembles the results obtained from the vertex method and significantly reduces overestimation resulting from calculations involving interval numbers.
In addition, it is worth noting that all calculations were conducted using MATLAB with the IntLab toolbox (Rump [60]), which provided essential support for interval arithmetic computations.

Conclusions
This paper presents a method for determining the dynamic response of systems with viscoelastic damping elements, where parameters are uncertain.It is assumed that uncertain parameters are described using interval numbers, meaning that the lower and upper bounds of the parameter are known.The presented method enables the consideration of uncertainties related to both the construction parameters and the damping elements, with a focus on demonstrating its applicability to parameters of viscoelastic material.
To determine the system's dynamic response, interval analysis was applied.However, direct application can lead to the significant overestimation of results.To reduce this overestimation, equations of motion were formulated using the element-by-element technique, and a transformation was applied, ensuring uncertainties related to the same parameter are considered only once.This dual approach effectively reduces overestimation, as demonstrated in the provided examples.Additionally, all results were compared with those from the vertex method, known for its accuracy but being computationally intensive.The method used here involves solving linear interval equations via the fixed-point iteration method, utilizing Brower's fixed-point theorem.It has been shown to yield satisfactory results at a considerably lower computational cost than the vertex method.However, it was observed that in cases of larger parameter variability (e.g., 10%), this iterative method may not always converge, necessitating the monitoring of convergence during computations.Nevertheless, as demonstrated Im(H44()) in the examples with variability of 5%, the method produces accurate results and can effectively estimate the lower and upper bounds of the dynamic response of structures.This method was applied for the first time to systems with viscoelastic elements described by both classical and fractional rheological models.Formulations were also proposed for both the Kelvin and Maxwell models.The Maxwell model, to be applicable with the presented method, requires a formulation with internal variables.The proposed approach is general and can be applied to more complex models, such as the Zener model or the generalized Maxwell model.
In future research, further reductions in overestimations, particularly in resonance zones, and improvements in method convergence in the case of larger uncertainties are planned.Additionally, the method is envisaged to be applied to the analysis of beams and plates with viscoelastic elements, which could have significant practical implications.
where      represents the sensitivity of the response function computed with respect to the changing parameter  , which lies on the boundary of the defined interval    ,   .The second term of the Taylor series represents the increment in the response function, which, after considering the lower and upper bounds of parameter  , can also be expressed as an interval number as follows: Δ  Δ  , Δ  (A7) Hence, the lower and upper bounds of the response function can be calculated as follows: (A8) To determine the lower and upper bounds of the increment for each changing parameter, it is necessary to determine the sensitivity of the response function and then select the smallest and largest increment values.This approach is detailed in [20], while the method for computing the sensitivity of the frequency response function can be found in [61].

Figure 3 .
Figure 3.A frame with two Kelvin dampers.
by Lewandowski et al.This formulation expands the dimension of the problem by the number of internal variables, yet it enables the

Table 1 .
Values of parameters for example 6.1.

Table 2 .
Ranges of parameters for the damper for example 6.1.

Table 3 .
Values of parameters for example 6.2.

Table 4 .
Ranges of parameters for the damper for example 6.2.

Table 5 .
Dynamic characteristics of the frame from example 6.3.

Table 9 .
Ranges of design parameter values.
(34)tional matrices  , ,  , ,  , ,  , ,  , , and  , are associated with additional degrees of freedom related to internal variables, while matrices  , and  , are associated with degrees of freedom of the analyzed structure (Lewandowski et al.[41]).The matrices are constructed using the element-by-element technique, and transformation(34)is applied.The results, in the form of selected elements of the FRF matrix, are shown in Figures24-27.

Table 10 .
Dynamic characteristics of the frame from example 6.4.