Identiﬁcation of Fractional Damping Parameters in Structural Dynamics Using Polynomial Chaos Expansion

: In order to analyze the dynamics of a structural problem accurately, a precise model of the structure, including an appropriate material description, is required. An important step within the modeling process is the correct determination of the model input parameters, e.g., loading conditions or material parameters. An accurate description of the damping characteristics is a complicated task, since many different effects have to be considered. An efﬁcient approach to model the material damping is the introduction of fractional derivatives in the constitutive relations of the material, since only a small number of parameters is required to represent the real damping behavior. In this paper, a novel method to determine the damping parameters of viscoelastic materials described by the so-called fractional Zener material model is proposed. The damping parameters are estimated by matching the Frequency Response Functions (FRF) of a virtual model, describing a beam-like structure, with experimental vibration data. Since this process is generally time-consuming, a surrogate modeling technique, named Polynomial Chaos Expansion (PCE), is combined with a semi-analytical computational technique, called the Numerical Assembly Technique (NAT), to reduce the computational cost. The presented approach is applied to an artiﬁcial material with well deﬁned parameters to show the accuracy and efﬁciency of the method. Additionally, vibration measurements are used to estimate the damping parameters of an aluminium rotor with low material damping, which can also be described by the fractional damping model.


Introduction
In modern engineering applications, advanced materials and material combinations are used to tackle vibration problems. In order to analyze such structures with a virtual model, an accurate representation of the structure, including geometric details, material characteristics, loading conditions etc., is required. While the geometric description is generally available from CAD data, the material parameters of the model are often unknown and have to be determined by comparison with measurement results.
For an appropriate description of the real structure, viscoelastic material behavior is often required. While a perfectly elastic material has three assumptions, namely the linearity, the simultaneity, and the unique equilibrium value between the stress and strain [1], viscoelastic material models only require linearity [1]. Rheological models for these materials can be illustrated as combinations of springs and dampers [2]. In the present paper, the fractional Zener model is used, which is a modification of the traditional Zener model presented in [3]. While the classical Zener model contains two springs and a dashpot [4], the introduction of fractional calculus leads to a modification of the dashpot behavior, see, e.g., [4,5], and an additional fourth parameter, the fractional-order of the time derivative, is introduced [6]. This modification allows for the representation of real materials with a low number of parameters [6].
A possibility to estimate the parameters in the material model is through a minimization of the error between a numerical model and measurements. Different experimental setups are available to generate reference data for the minimization process. Quasi-static methods, such as the creep test and the relaxation test, analyze the behavior of the materials to an impulsive loading over a long time [1], while dynamic methods consider the response to varying loads. According to [7], the dynamic methods can further be categorized into the Dynamic Mechanical Analysis (DMA) [8] and vibration tests [9]. Several other techniques are listed in [9], with a special focus on the estimation of damping.
The analysis of the material behavior by DMA is based on the measurement of the amplitude and the lag in the stress-strain relation caused by an oscillating force [8]. Then, the parameters of the model, such as the modulus and the damping, can be calculated [8].
In [10], the viscoelastic material parameters of synthetic rubber are estimated by DMA. The damping of the test object is described by a Maxwell model, where two different specifications are compared [10]. First, a series Maxwell setup containing 12 components and second, a fractional derivative model containing a springpot element instead of the dashpot is analyzed [10]. In [11], a unidirectional glass fibre-reinforced epoxy is investigated, where a three-point bending test including the effect of different temperatures of the test object is applied.
DMA leads to very good estimates for the material parameters, but expensive test equipment is required [12] and the method is limited to low frequencies [7]. Alternatively, vibration tests can be used to estimate the material properties. In [7], a bare beam and a sandwich structure are analyzed by shaker tests. The Frequency Response Function (FRF) curve is used in a model updating process applying the amplitude correlation coefficient [7]. The results are compared to those of DMA and show that the vibration test was found superior [7]. In [12], a setup containing two beams with different lengths is analyzed by impact hammer tests. The response is measured by an acceleration sensor and the complex Young's modulus is estimated, based on the comparison of a calculated FRF of the beam and the measured response [12].
Comparing the measured FRF curves with numerical data is also used by other researchers, like in [13][14][15]. In [13], isotropic and orthotropic material parameters are estimated from bending tests. Finite element (FE) simulations predict the viscoelastic material behavior and the Levenberg-Marquard algorithm is used for the estimation of the best fitting parameters [13]. In [14], a finite element model is applied to describe the test object containing a beam with a layer of viscoelastic damping material in a clamped-free boundary condition. The FRF curves are used to extract the values of the eigenfrequencies and amplitudes at these eigenfrequencies for the Parameter Identification Process [14]. A two-step minimization is applied to estimate the best fitting parameters in combination with a parameter sensitivity analysis [14]. In [15], a viscoelastic composite plate is investigated. Amplitude values at the resonance frequency and close to the resonance frequency are used to fit the parameters [15].
Generally, simple objects such as beams or plates are used in the process of material identification, since the modeling is less involved and the measurements are easier to conduct. Nowadays, various methods exist to solve one-dimensional structural problems, such as the Transfer Matrix Method, the Dynamic Stiffness Method or the Green function method [16]. These methods involve a frequency-dependent system matrix that needs to be solved numerically [16]. In the present paper, a (semi-)analytical method, named Numerical Assembly Technique (NAT), is used to describe one-dimensional beam structures. The method was introduced by Wu and Chou [17] in 1999. NAT subdivides the beam structure into segments that contain a constant cross-section and constant material parameters [16]. The homogeneous governing equations of the uniform beam segments are solved analytically and the resulting solutions are used to fit the boundary and interface conditions [16].
Generally, the minimization process for the parameter identification requires numerous evaluations of the numerical models, which is a time-consuming process. In order to raise the efficiency of the Parameter Identification Process, the numerical model can be replaced by a surrogate model. Several surrogate modeling techniques exist, such as the Artificial Neural Network [18][19][20][21], parametric Reduced Order Model [22], the Response Surface Method [23] or the Polynomial Chaos Expansion (PCE) [24][25][26][27][28]. In [18,19], a viscoelastic cylinder coated by water is analyzed and the complex Young's modulus is estimated. An Artificial Neural Network is used as a surrogate model. In [20,21], a laminated structure is analyzed. The FRF curve gives the reference values and an Artificial Neural Network is used for the surrogate modeling process [20,21]. A Bayesian approach is applied for the estimation process [20,21]. In [22], a reduction in the computational time is achieved by the reduction in the order of the model in order to efficiently estimate the complex shear modulus. In [23], a Response Surface Method is used to describe the reference values based on the storage modulus and the loss factor. In [24], a frictional system is analyzed and the Polynomial Chaos Expansion is used as a surrogate modeling technique. In [25,26], composite plates are analyzed and, based on the deviation of the results, the input parameters are chosen to describe the polynomial basis. In [27,28], a special combination of the Polynomial Chaos Expansion to estimate the material parameters [27] and bearing stiffness parameters of a test rig [28] is shown. The used method is called Polynomial Chaos Kriging and uses the Polynomial Chaos Expansion for a global system description and the Kriging method for the interpolation of local values [29].
In the present work, PCE is used, which approximates the complex numerical model by orthonormal polynomials to reduce the computational load. PCE is based on the concept of Wiener [30] and has been extended by Cameron and Martin [31] using nonlinear functions defined by a Fourier-Hermite series [32][33][34]. Many refinements and new fields of applications for PCE have been investigated, which are summarized by Ghanem and Spanos [32]. The aim of PCE is to span a space based on a polynomial basis and correlating coefficients in order to build a surrogate model [34]. Depending on the distribution of the input parameters, different polynomials can be applied, as listed in, e.g., [33].
The outline of the paper reads as follows: In Section 2, the two numerical methods are explained, which build the basis for the Parameter Identification Process (PIP). In Section 3, two examples are presented. First, a numerical example shows the efficiency and simplicity of the presented method. Second, a real measurement is used and the material and damping parameters are estimated. Finally, in Section 4, the conclusion is presented.

Materials and Methods
In this section, the Numerical Assembly Technique (NAT) and the Polynomial Chaos Expansion (PCE) are explained. These methods form the basis of the Parameter Identification Process (PIP).
NAT is an efficient computational technique, which is used to analyze the dynamic behavior of beams. A numerical model includes input parameters, such as the material parameters and geometry parameters. Based on these input parameters and an excitation function, the dynamic response of the system is calculated. This gives the FRF curve, from which parameters such as the eigenfrequency of the beam are extracted. These parameters are named output parameters in the following. So, NAT describes the correlation between input parameters and output parameters.
The aim of the PIP is to find appropriate values of the input parameters for given output parameters. In order to find these input values, the model needs to be evaluated very often, which is very time-consuming. Here, the computational effort is reduced by a surrogate model: PCE uses polynomials to describe the correlation of the input parameters and output parameters to raise the efficiency of the estimation process.

Numerical Assembly Technique (NAT)
The Numerical Assembly Technique gives a (semi)-analytical solution of one-dimensional structures [16]. In this section, NAT is explained for the Timoshenko beam theory with viscoelastic material behavior represented by the fractional Zener model. The content and the structure of the description is based on the work of the authors presented in [35].
NAT is based on the subdivision of a one-dimensional structure into segments. Each segment has a constant cross-section area and constant material parameters. The point between two segments is named station. At the stations, lumped masses, mass moments of inertia, springs or dampers can be applied. The governing equation describing a segment is solved. The resulting displacement, rotation, shear force and moment are parameters, which are used to fulfil the boundary and interface conditions of the segment. All segments are assembled to give a full description of the structure [16,35].
The harmonic equations of the Timoshenko beam theory are given by [35] where j is the imaginary number, the index is the segment number,w (x) is the complex amplitude of the transverse displacement w (x, t) = Re[w (x) e jωt ],φ (x) is the complex amplitude of the rotation ϕ (x, t) = Re[φ (x) e jωt ],M (x) is the complex amplitude of the moment M (x, t) = Re[M (x) e jωt ],Q (x) is the complex amplitude of the shear force Q (x, t) = Re[Q (x)e jωt ], x is the global coordinate, ω is the angular frequency, ρ is the density, A is the cross-section, I is the second moment of area, k S is the shear correction factor, d a is a damping coefficient,q(x) is the complex amplitude of the distributed is the real part of •, E is the complex Young's modulus and G is the complex shear modulus. The complex Young's modulus E (ω) and shear modulus G (ω) are described by a four-parameter model called fractional Zener model. This model describes a viscoelastic material behavior. The four-parameter model was analyzed by Pritz [6]. The axial stress based on this model is described by [35] whereσ xx (x, z, ω) is the normal stress,ε xx (x, z, w) is the normal strain and a E 0 , and a E 1 , b E 0 and α E are four positive real constants. Further, the shear stress is given by [35] whereσ xz (x, ω) is the shear stress,γ xz (x, z, w) is the shear strain, and a G 0 , a G 1 , b G 0 and α G are four positive real constants [35].
The numerical description of the system given byw (x ) is represented by a homogeneous solutionw h (x ) and a particular solutionw p (x ) of the governing equation. First, Equation (1) is solved by neglecting the excitation. Considering the boundary conditions given in [36] and the modification ofw h = c w e jkx , wherew h is the homogeneous part of the displacement, a system of equations representing the homogeneous solution is given by [35]x Here,x h (x ) contains the homogeneous field variables, B (x ) is a matrix containing the material parameters and c is a vector of constants. Second, the particular solution is calculated, which considers the effect of the excitation on the particular solution ofx p (x ) [16].
Last, the assembly process is operated where a system is described by [16] where A is a system matrix containing the matrices B and the vectorc accounts for the constant parameter values which are to be estimated. The vector b contains the loadings and boundary conditions [16].

Polynomial Chaos Expansion (PCE)
The Polynomial Chaos Expansion is a surrogate modeling technique which has many applications. In the present case, the method gives a very efficient surrogate model of the structure. More investigations on this method are given in [34,37].
The content and the structure of the following are similar as given in [37]. A quantity of interest Y with a finite variance is described by an infinite summation given by [37,38] where y i represents coefficients and Z i represents a numerable set of random variables. Instead of the numerable set of random variables, PCE uses a set of multivariate orthonormal polynomials. For this, a polynomial basis needs to be defined, which represents the multivariate orthonormal polynomials. The construction of the basis starts with univariate orthogonal polynomials π k represents the orthogonal polynomial, k represents the degree of the polynomial and ζ i represents the independent variable. The orthogonality is described by [37] where E[] represents the expectation value, δ jk is the Kronecker delta and a i j is the squared norm of the polynomials given by [37] This orthogonal polynomial basis is usually not orthonormal. Therefore, the orthogonal polynomials are normalized by [37] where ψ (i) j is the orthonormal polynomial. A table of orthonormalized polynomials is given in [37]. The polynomial basis of the Polynomial Chaos Expansion is built up by the univariate orthonormal polynomials by [37] where ψ α (ζ) represents the multivariate orthonormal polynomial and α the multi-indices.
The multi-indices are ordered lists of integral. Finally, the Polynomial Chaos Expansion is described by [37] where Y represents the surrogate model and y are the corresponding coefficients of the polynomials ψ α which need to be estimated in order to describe the system. In practice, the infinite series expansion is truncated. The truncated series is given by [37,39] where A is the truncation factor. This factor is defined by [37] card where M is the number of independent input variables, p is the polynomial order and •! is the factorial function [37]. Given the description of the multivariate orthonormal polynomial basis, the coefficients need to be estimated. Two categories of computational schemes were developed to estimate the coefficients: the intrusive schemes and the non-intrusive schemes [37].
Intrusive schemes are applied, e.g., in the stochastic finite element method. Equations, which describe the system, are discretized in the physical space and in the random space. The results are coupled and solved intrusive. The non-intrusive schemes use a realization of the model by repeated runs and a, e.g., least-square minimization in order to estimate the coefficients [37].
The aim of this work is to describe the material constants given in Equations (6) and (7) in order to achieve this in an efficient way, the system as given in Equation (9) is represented by a surrogate model as given by Equation (16). This constitutes the model for the Parameter Identification Process.

Parameter Identification Process (PIP)
The structural model presented in Section 2.1 contains several input parameters. These parameters describe the geometry of the structure and the material properties. Some of these parameters are easy to measure, such as the length of a section or the density of the material. In the following, a new process is developed to estimate the parameters of the fractional derivative model. This Parameter Identification Process (PIP) uses reference values from a measurement, which have to be met by the structural model using the estimated input parameters.

Reference Values and Fundamentals of the PIP
The measured FRF is used to estimate the reference values of the identification process. Peaks in the FRF curve are determined, where the eigenfrequency (EF), the amplitude at the eigenfrequency (AMP) and the normalized frequency band (NFB) are extracted as the reference values for each peak. The NFB is equivalent to the amplification factor given in [40]. This factor is estimated by the half-power method [40]. In [41], the NFB value is called the quality factor of a system and its reciprocal value is also named loss angle [1], which is an indicator for the internal friction of the material. The NFB value is given by [1,40] The frequencies f 1 and f 2 are defined by the crossing point of the Frequency Response Function curve and a horizontal line which marks 70.7% of the amplitude at the eigenfrequency. The 70.7% value is often used for the description of the structural or hysteretic damping of the system and can be found in different sources, like in [40] or [1]. In Figure 1, a peak of an FRF curve including the relevant reference values is illustrated. The presented PIP is shown for homogeneous isotropic materials. At the beginning, the complex Young's modulus E (ω) is rewritten as [1,42] where E (ω) is the frequency depending storage Young's modulus and E (ω) is the frequency depending loss Young's modulus [1,42]. These two parameters are combined to |E (ω)| given by where |E (ω)| is the frequency-dependent absolute value of the complex Young's modulus. The parameter η(ω) is given as the ratio of the loss and the storage modulus [1,42] The value η(ω) describes the frequency depending energy loss of the material [1,42]. This allows for rewriting the complex Young's modulus as The reference values estimated from the FRF curve are values depending on a specific frequency. The complex Young's modulus as represented in Equation (22) describes the value over a frequency range. For the following process, the unknown parameters are estimated at the specific frequencies defined by the frequency of the reference values.

Bisection Method
The PIP focuses first on the estimation of |E (ω)| and η(ω). Other parameters related to the mass or the geometry of the structure are assumed to be known or easy to measure. The first step of the process uses bisectioning on a wide range of parameter values. This method is basically used to estimate the roots of a function and uses the halving of the analyzed interval [43]. The limit values of |E | and η and the mean values are defined. This gives 9 possible parameter combinations for the bisectioning as depicted in Figure 2. Only four out of these combinations are evaluated: For the combinations marked in grey color, the eigenfrequencies and corresponding amplitudes are computed from the NATmodel described in Section 2.1 and compared to the corresponding measured values at the specific frequency. The comparison leads to different cases to be considered: • The combination of |E µ | with both, η min and η max gives computed eigenfrequencies higher than the measured value. In this case |E max | is disregarded and replaced by |E µ |; • Similarly, if both combinations of |E µ | give a value of the computed eigenfrequency lower than the measured one, then |E min | is replaced by |E µ |; • For the combinations of η µ with both |E min | and |E max |, the computed amplitudes are compared to the measured amplitude. If both computed values are higher than the measured values, then η min is replaced by η µ ; • Similarly, if both combinations are lower than the measured amplitude, η max is replaced by η µ .
The bisection method stops, when defined limit values are reached or if one calculated value is higher than the correlating reference value and one calculated value is lower than the correlating reference value. The usage of the eigenfrequency and the amplitude as reference values is adapted from [14], where the estimation of the unknown parameters is explained as the two-step identification method. In [14], the reference values EF and AMP are proved successfully for the fitting process. The third reference value is used in the bisection method as an indicator of the correct estimation of η, because [1] The bisectioning gives rough estimates for |E | and η, which are used to generate a surrogate model.

Surrogate Modeling
The surrogate model describes the relationship between the input parameters |E | and η on the one hand side and the output values EF, AMP and NFB on the other side.
For each output parameter, an independent surrogate model is generated. Depending on the number N of peaks in the FRF curve, 3N surrogate models are built. In a general case, no information is given about the distribution of the input parameters, so the input parameters are chosen to be uniformly distributed. After determining the surrogate models, a large input sampling set is generated for each peak in the FRF curve depending on its limit values.
The surrogate models are evaluated for the values given by the sampling set and the corresponding output parameters are calculated. The differences between calculated and measured values are analyzed and based on this information, new limit values (the final limit values) for the input parameters |E | and η are defined. With these final limit values, new surrogate models are generated for each reference value. These final surrogate models are estimated and the corresponding final limit values for each peak in the FRF curve are determined.

Curve Fit
After the estimation of the final surrogate models and its corresponding limit values of |E | and η, a large input sample set is generated and the surrogate models calculate the output parameters EF, AMP and NFB for each peak in the FRF curve. These calculated values per peak are compared to the reference values per peak and a finite number of fitting parameter combinations of |E | and η are estimated. These finite number of parameters are converted into the real and imaginary parts of the complex Young's modulus E . This conversion has the advantage that two curves are available for the non-linear curve fit for the description of the complex Young's modulus over the frequency range. The curve fit is executed based on the number of supporting points, which is equal to the number of peaks in the FRF curve. In Figure 3a,b the real and imaginary parts of the complex Young's modulus for an numerical example are illustrated.  The complex Young's modulus described in Equation (22) is describable with fractional derivative parameters based on the description of a fractional Zener model given by [35] Following, the support points for the description of the real and imaginary parts of the complex Young's modulus are used to fit the fractional derivative parameters presented in Equation (24).

Estimation of the Fractional Derivative Parameters
After the curve fitting process, possible fractional derivative parameter combinations are estimated. Based on the maximum error of the curve fit of the imaginary and the real part of the complex Young's modulus, a set of appropriate fractional derivative parameter combinations is estimated. Based on these limit values of the input parameters, a large sample set of input parameters is generated and the output parameters EF, AMP and NFB are calculated from the NAT model. The values that fit the measurement are used for the generation of the surrogate models. After the generation of the surrogate models, the best fitting fractional derivative parameters are estimated.
Therefore, the range of the values of the possible fractional derivative parameters is equidistantly subdivided. Under the consideration of [35] a large number of input sample sets is generated and the surrogate models calculate the output parameters. Following, the error of the input sample set is estimated and an error minimization is given by where REF represents the reference values EF, AMP and NFB. Based on the minimum error, the global optimum is estimated by refining the input sample space and the repeated calculation of the output parameters based on the surrogate model. Closing, the best fitting combinations for the surrogate model are estimated. These combinations are solved with the numerical model and the FRF curve of the measurement and the FRF curve based on the fractional derivative parameters estimated by the PIP are compared.

Results
In this section, the PIP technique is applied to estimate the fractional derivative parameters of an unknown material. First, an example is illustrated step by step to show the efficiency and accuracy of the method. Therefore, a theoretical material is used to generate a FRF curve and following, the fractional derivative parameters are estimated. Second, the process is applied to a real structure. A beam made of aluminium is analyzed and the fractional derivative parameters are estimated. NAT is implemented in MATLAB ® 2020b and the surrogate modeling process is executed by the MATLAB ® toolbox UQLab [44].

Numerical Experiment
A numerical example is used to show the PIP in detail. The beam is made of a nonexisting material described by the fractional Zener model. These reference parameters are listed in Table 1. First, the geometry of the system needs to be defined. Here, a beam with six stations and five segments is considered. The cross-sections of the beam are modeled as being circular. The diameter d 1 = d 2 = d 5 = 0.04 m and the diameter d 3 = d 4 = 0.05 m describe the beam. The density of the beam is ρ = 8250 kg m 3 and the Poisson number is ν = 0.27. The shear correction factor is calculated based on [45] for a homogeneous circular cross-section as [45] k S = 6 · (1 + ν) 2 7 + 14ν + 8ν 2 = 0.8516, (27) where k S is the shear correction factor. The position of the stations, the additional mass moments of inertia Θ, and additional masses m are listed in Table 2. In Figure 4, the numerical example is illustrated. The FRF curve is calculated in the case of the numerical example, and all the reference values are extracted. Figure 5 illustrates the FRF curve that represents a point force excitation at x F = 0.34 m and a response at x R = 0.015 m.  As a first step, the peaks in the FRF curve (represented by dots in Figure 5), containing the reference values EF, AMP and NFB need to be determined. All four peaks are used to describe the fractional derivative parameters of the system. In Table 3, the reference values are listed. Based on the resolution of 1 Hz, the frequency values needed for the estimation of the NFB value are calculated by a linear interpolation. After the estimation of the reference values, the limit values of |E | and η of each peak in the FRF curve are determined from the bisection method presented in Section 2.3.2.
The limit values of |E | and η are based on the EF and AMP reference values. The bisectioning gives rough limit values, where a difference of 10% to the corresponding EF reference value and 20% to the corresponding AMP reference value per peak in the FRF curve is acceptable. These tolerances are a rough limit and can be defined freely, e.g., with reference to the measurement equipment. In Table 4, the limit values of |E •,B | and η •,B depending on the peak in the FRF curve after the bisection processes are listed. By the calculation of the reciprocal value of the NFB value, the correctness of the limits of η •,B is shown, referred to Equation (23). Based on this check, the correctness of the estimated values of η •,B is proofed and the PIP is continued.
In the next step, the surrogate models are generated with the MATLAB ® toolbox UQLab [44] based on the values listed in Table 4. The estimated limit values are used as the limit values of the uniformly distributed input parameter. Legendre Polynomials are applied to generate the multivariate orthonormal space [37]. The polynomial order p is defined and the error between the calculated results by NAT and the estimated values by the surrogate model is minimized. Based on the polynomial order p and the number of independent variables M, Equation (17) calculates the number of needed output parameters for the estimation of the surrogate model.
These generated surrogate models are evaluated and based on the quality of the model, new limit values of |E •,S | and η •,S are defined. These new and final limit values are listed in Table 5. Again, the reciprocal values of the NFB values listed in Table 3 are between in the limit values of the η •,S values. This shows the correct estimation of the parameters. The final surrogate models, as described in Section 2.3.3, are used for the estimation of the possible parameter combinations for the curve fitting process. Large sample sets, based on the limit values listed in Table 5 are generated and a finite number of possible parameter combinations per peak in the FRF curve are estimated. These finite number of parameter combinations are now transformed into the real and the imaginary part of the complex Young's modulus. Closing, the non-linear curve fit is executed for the estimation of the fractional derivative parameters described by Equation (24).
The non-linear curve fit process is executed with the MATLAB ® function lsqcurvefit. Based on these non-linear curve fits, the parameter combinations of the fractional derivative, based on the description of the complex Young's modulus, is given. Afterwards, the normalized error at each peak in the FRF curve related to its reference value of the real part and the imaginary part is estimated.
In Figure 6a,b the errors normalized to the real and imaginary parts of the parameters estimated from the surrogate modeling are illustrated. Based on the lowest errors, the possible fractional derivative parameter values are estimated. In Table 6  The presented surrogate modeling process and error minimization process gives an estimation for the best fitting parameter combination. Table 7 lists the PIP-approximated parameters.  Figure 7, the measured and calculated FRF curves of the numerical example are illustrated. In [6], the quality of a the parameter fit is evaluated by the comparison of the static Young's modulus a E 0 . The value a E 0,Re f = 7 × 10 10 N m 2 is used to generate this example; the value a E 0,PIP = 7.0036 × 10 10 N m 2 is found by PIP. The relativ error is 0.051%. This shows the efficiency of the PIP. Additional, η(ω) is analyzed, to illustrate the difference between the parameters that are used for the generation of this example and the PIPestimated values.
In Figure 8, the frequency-dependent loss value η(ω) of the given values of the numerical example and the fractional derivative parameters estimated by the PIP are illustrated. The dots represent the frequency points used for the PIP, which are illustrated by the peaks in the FRF illustrated in Figure 5. It is shown that with four peaks as reference, the system is described accurately.  To assess the efficient of the PIP, the computational time to compute the training set, to generate the surrogate models and the actual Parameter Identification Process is compared to a straight forward identification process without a surrogate model. All calculations are carried out on a computer operating on Windows 10, with an Intel ® Xeon ® E3-1270 processor (4 × 3.6 GHz) with 32 GB RAM. The calculations are executed with MATLAB ® 2020b.
The first step, the bisectioning process, only requires a small number of calculated reference values; therefore, the computational time is negligible compared to the other steps. The surrogate modeling process uses 8000 reference calculations to estimate the limit values of |E | and η at each reference frequency. The total computational time to calculate these reference values is approximately 6800 s. Depending on the used polynomial order, the generation of the surrogate model takes between 0.07 s (polynomial order p = 3) and 2.2 s (polynomial order of p = 30). The needed parameter combinations for the description of the metamodel are calculated by Equation (17), the remaining parameter combinations are used for the quality check of the PCE. As an example, for the computation of 10 6 input parameter sets (required to find a reasonable estimate of |E | and η), the original NAT model requires about 236 h, while the total calculation time of the surrogate model including the training process is only 6826 s. Therefore, the surrogate model approach is 124 times faster in this case.

Measurements
In the following, the PIP is applied for measured FRF curves obtained from a real structure.
The analyzed beam is made of aluminium. During the measurement, the test object is suspended horizontally by fishing lines with a diameter of 1 mm, giving a free-free boundary condition. The test setup is illustrated in Figure 9.
The beam is excited with an electrodynamic shaker (Brüel and Kjaer LDS V406, Virum, Denmark) and the response is measured with a triaxial acceleration sensor (Brüel and Kjaer 84506, Virum, Denmark). The excitation force is measured with a force transducer (Brüel and Kjaer 8230-001, Virum, Denmark), where a stinger is mounted between the shaker and the force sensor. The sensors are installed on the test object with a thin layer of wax on the surface. Both sensors are in the same horizontal plane. The excitation signal is sinusoidal and the system is in steady state condition when the response is measured. The duration of the measurement is one second; the resolution of the measurement is 1 Hz. First, the time signals are recorded. Then, a FFT is obtained from the measured data and the dynamic response of the test object is estimated. Finally the Frequency Response Function signal is transformed into dB by the transformation of the logarithmic function. In Figure 9, the multi-stepped beam made of aluminium a is illustrated. The force excitation is applied at the global position x F = 340 mm, in which the force is measured by a force transducer b . This sensor is connected with a stinger c to a electrodynamic shaker d . The response is measured at the global position x R = 565 mm by a triaxial acceleration sensor e . The support, realized by fishing lines f , represents the free-free boundary condition. Screws g are used to fix the fishing lines at the beam. In Figure 10, the geometry of test object, including the position of the lumped masses and mass moments of inertia, is presented.
The stepped beam is subdivided into five segments, in which the cross-section is described by a circular and defined by diameters d 1 = d 4 = d 5 = 0.04 m and d 2 = d 3 = 0.05 m. The density is calculated based on the measured mass and the calculated volume of the beam as is given by ρ = 2798 kg m 3 . The Poisson ratio is based on a literature value of ν = 0.34 [46]. Equation (27) gives the shear correction factor based on the Poisson ratio, which is given by k S = 0.8493. Table 8 lists the position of the stations, the lumped masses and the mass moment of inertia of the system.  The additional masses and mass moments of inertia represent the acceleration sensor (station 5) and the force transducer (station 3) of the system, in which the lumped mass (station 1 and station 6) represent the mass of the screw.
The measured FRF curves are illustrated in Figure 11. Four curves are measured, while the test configuration is identical. The modification between the measurements is based on the variation of the force amplifier to analyze any modifications. It is shown that the measurement values are nearly equal. The anti-resonance at around 1000 Hz illustrates the different amplifier levels of the measurements.  The reference values of the measurement define the limit values of the PIP. Based on four different measurements, a range of the reference values per peak in the FRF curve is given. The limit values are listed in Table 9. Based on the resolution of the measured Frequency Response Function curve of 1 Hz, the values that are needed for the estimation of the NFB are linearly interpolated. It is shown that the variation of the eigenfrequencies EF and the amplitude AMP is small, whereas the variation of the Normalized Frequency Band NFB is larger. This is caused by the dependency of the NFB by the AMP. Applying the bisection part, the surrogate modeling part and the curve fitting of the PIP, the limit values for estimating the fractional derivative parameters are determined.
The limit values are listed in Table 10. Finally, the best fitting parameters are determined under the consideration of global minimization. Based on the range of the limit values, a range of fitting parameters is estimated. One possible parameter combination is listed in Table 11. Therefore, one measurement is used for global minimization.   The quality of the modeling process is shown by the accuracy of the anti-resonance around 1000 Hz. The value of a E 0 is similar to the value listed in the literature [46] for the static modulus of aluminium. It is shown that the system is fitted nearly perfectly. The maximum error is 1%. This error is shown at the third bending frequency.

Conclusions
The presented method shows the efficiency of combining a one-dimensional analytical description of structures with high-order polynomial surrogate modeling. No information about the system is needed initially and the limit values are estimated fast and efficiently. Based on the NFB value, a quality check of the estimated values is executed. Using surrogate modeling of the fractional derivative parameters, an estimate of the parameters can be found which describes the whole structure over a wide frequency range.
The advantages of the presented method are as follows: • Analytical solutions of real structures are indispensable for a minimization of the error between the mathematical description of viscoelastic behavior and the real behavior of materials. • The Polynomial Chaos Expansion, especially the MATLAB ® toolbox UQLab [44], enables a fast and efficient implementation of surrogate modeling. The high-order polynomials are an efficient way to surrogate a large numerical calculation with a minimized error. • Based on splitting up the complex Young's modulus into the real part and the imaginary part, two curves are used for the parameter fit, which minimizes the number of needed peaks in the FRF curve. • The final surrogate modeling process is used to describe the whole structure with a minimum error in the representation of the whole system and the equidistant split of the input sample rang the global minimum is found directly. Funding: Open Access Funding by the Graz University of Technology.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author. The data are not publicly available due to an ongoing research project.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript: