Steady-State Harmonic Vibrations of Viscoelastic Timoshenko Beams with Fractional Derivative Damping Models

: Due to growing demands on newly developed products concerning their weight, sound emission, etc., advanced materials are introduced in the product designs. The modeling of these materials is an important task, and a very promising approach to capture the viscoelastic behavior of a broad class of materials are fractional time derivative operators, since only a small number of parameters is required to ﬁt measurement data. The fractional differential operator in the constitutive equations introduces additional challenges in the solution process of structural models, e.g., beams or plates. Therefore, a highly efﬁcient computational method called Numerical Assembly Technique is proposed in this paper to tackle general beam vibration problems governed by the Timoshenko beam theory and the fractional Zener material model. A general framework is presented, which allows for the modeling of multi-span beams with general linear supports, rigid attachments, and arbitrarily distributed force and moment loading. The efﬁciency and accuracy of the method is shown in comparison to the Finite Element Method. Additionally, a validation with experimental results for beam systems made of steel and polyvinyl chloride is presented, to illustrate the advantages of the proposed method and the material model.


Introduction
In the modeling process of real structures, many components are accurately described by simplified beam models, e.g., Euler-Bernoulli beam or Timoshenko beam theory, if two dimensions are significantly smaller than the third one. Therefore, a detailed analysis of this simplified structural elements is very important for practical applications; hence, a vast amount of literature on this topic is available, and many efforts have been devoted to understand the dynamics of beam models [1]. Due to advanced manufacturing processes, more complicated materials, especially viscoelastic materials [2], can be used in the product design to further improve dynamic characteristics, e.g., noise emission or vibration levels [3]. Therefore, it is necessary to also include such materials in the analysis of beam vibrations.
An accurate modeling of the behavior of such materials is required and various constitutive equations have been presented, e.g., the Kelvin-Voigt model [4], the Maxwell model [5], or the Zener model [6]. Increased possibilities in mechanical testing have shown that many materials with complex microstructure lead to a power-law signature in their creep and relaxation behaviors [7]. Even though recent studies, e.g., References [8,9], show that generalized versions of the classical models can lead to a good agreement with measured data for such materials, a large number of model parameters is required, which greatly hinders physical interpretation [7]. Therefore, so-called fractional time derivatives have been introduced in the constitutive equations, which are ideally suitable to model the viscoelastic behavior of materials [7,10] and allow for an accurate representation of the material behavior with a small number of parameters over a wide frequency range [11,12]. particular solutions are derived. A numerical validation of the proposed method is given in Section 4 using a numerical reference solution computed with the Finite Element Method. Additionally, a comparison with vibration measurements on two different beam structures made of steel and PVC is carried out. Finally, a conclusion is given in Section 5.

Problem Description
In this section, a general two-dimensional beam vibration problem is described, and different viscoelastic material models using fractional time derivatives are outlined. The Timoshenko beam theory is applied to model the beam segments including external viscous (air) damping.

General Viscoelastic Beam Vibration Problem
A general beam vibration problem is illustrated in Figure 1. It is assumed that the system is two-dimensional; therefore, no deformation occurs in the y-direction. Additionally, the system has no loading in the x-direction. The beam with total length L is split into M uniform and homogenous segments by (N) stations. The first (1) and last station (N) represent the boundaries of the beam system. An additional station (i) has to be introduced if a discontinuity in the system appears, either through a concentrated element, a change in the cross-section or the material parameters. The global coordinate of an intermediate station (i) is given by x = X i and for each segment , a local coordinate system (O , x , y , z ) is defined. The origins of the local coordinate systems O are coincident with the bending centers of each beam segment and are located at the associated left station (x = (x − X i ) (i = )). Each beam segment is straight and uniform with the length L = (X i+1 − X i ) (i = ), complex Young's modulus E , complex shear modulus G , density ρ , cross-section area A , and second moment of area about the y-axis I . It is assumed that the principal axes of all segments are aligned; therefore, the axial and bending deformations are decoupled [34].
The support at station (i) of the beam is modeled by linear springs and dampers (translational springs k r ) and rigid attachments by lumped masses m (i) , and rotatory inertia Θ (i) . The translational and rotational springs and dampers are linear, and, in the initial undeformed state, the springs are unstressed. All concentrated elements are connected to the beam at the bending center line. As illustrated in Figure 1, it is assumed that a general support, a rigid attachment, and a change in material and geometric properties appear simultaneously at each station. Every possible configuration at a station (i) can be achieved by setting certain parameters of the supports or rigid attachments to zero or by keeping the material or geometric properties unchanged.
The external loading of the beam is modeled by point forces F P (t) at the local positions X F , point moments M P (t) at the local positions X M , distributed forces q(x, t), and distributed moments m(x, t). For steady-state harmonic vibrations at angular frequency ω, the external loads can be defined by where t is the time, j = √ −1 the imaginary unit, andF P ,M P ,q(x), andm(x) are the complex amplitudes of the concentrated and distributed external loads. The external viscous (air) damping is denoted by q D (x, t).

Modeling of the Internal Material Damping
Although many solid materials are modeled as perfectly elastic, mainly by Hooke's law, actual solids behave differently even for small stresses [35]. There is always at least a small dissipation of energy, even for metals [35].
In most beam theories, it is assumed that the transverse stress components σ yy and σ zz are very small compared to the axial stress component σ xx (σ xx σ yy , σ xx σ zz ). Consequently, it is assumed that the transverse stress components vanish (σ yy ≈ 0, σ zz ≈ 0) [34]. Furthermore, only small deformations are considered; therefore, a linear material model can be applied. A perfectly elastic material is defined by Hooke's law, which leads to the non-zero normal and shear stresses in a beam segment [34]: with E the Young's modulus, G the shear modulus, ε xx (x, z, t) the normal strain and γ xz (x, t) the shear strain. A transformation of Equations (2) and (3) into the frequency domain applying the Fourier transform results iñ A common approach to model the dissipative behavior of solid materials in the frequency domain is the constant hysteric damping model [36]. In this frequency independent model, the real-valued Young's and shear modulus are simply replaced by complex moduli, which leads toσ with the complex Young's modulus E = E (1 + j η E ) and the complex shear modulus [36]. The loss factor η • defines the frequency-independent damping in this model. Although this simple model has certain advantages for numerical simulation, especially in finite element schemes, and approximates experimental results accurately in a limited frequency band, it leads to non-causal responses in the time-domain [36]. Therefore, more realistic damping models have been developed. Several of these models are based on so-called fractional time derivatives, which can be used in a wide frequency range and for a large class of materials, such as metals, glass, and especially polymers [11]. A very effective four-parameter model has been analyzed by Pritz [11], which is a generalized form of the so-called Zener material model. The constitutive equations of this model in the time domain are given by [11] with a E 0 , a E 1 , b E 0 , a G 0 , a G 1 , and b G 0 as positive real constants, and the restrictions a E 0 < ∂ t α as the fractional derivative of order 0 < α < 1. Different mathematical definitions of the fractional derivatives have been presented in literature. In this paper, the Riemann-Liouville definition is applied, which has the property that the Fourier transform F {•} of the fractional derivative [37] allows for a simple representation of the four-parameter model in the frequency domaiñ The complex moduli E (ω) and G (ω) are frequency dependent. This model can be used to describe materials over a wide frequency range as long as the resulting loss factors η E (ω) and η G (ω) exhibit only one peak in the frequency range [11]. If several peaks are present in the loss factor, more evolved models have to be applied. The fractional Zener model is also capable of representing a low loss factor, which is practically frequency independent and, due to its simplicity and causality, superior to the hysteric damping model [11]. Therefore, the four-parameter model is used in the subsequent sections to model the material behavior of the beam segments . The four-parameter model can be reduced to simpler material models by setting certain parameters to specific values, e.g., if a • 1 and b • 0 is set to zero, Hooke's law is recovered, or b • 0 = 0 and α • = 1 leads to the Kelvin-Voigt material model.

Timoshenko Beam Theory
The beam segments are modeled by the well known Timoshenko beam theory, which can be found in several textbooks, e.g., References [34,38]. The displacements u (x, z, t), v (x, z, t), and w (x, t) of an arbitrary point C in x-, y-, and z-direction and the rotation of the cross-section ϕ (x, t) about the y-axis are shown in Figure 2a.
In Figure 2b, an infinitesimal element of a Timoshenko beam is shown. The equilibrium of forces and moments leads to with δ(•) the dirac delta function, the external viscous (air) damping with d a > 0 the damping coefficient, the bending moment and the shear force For harmonic vibrations at angular frequency ω, the state within a beam segment is completely described by the transverse displacementw (x, ω), the rotation of the cross-sectionφ (x, ω), the bending momentM (x, ω), and the shear forceQ (x, ω). Plugging Equations (11) and (12) into Equations (16) and (17) and, subsequently, into Equations (13) and (14) leads to and k S the shear correction factor to account for the actually shear stress distribution. For brevity, the dependency of the variables on ω in Equations (18)-(21) is dropped, and point loadings are not explicitly stated in the equations but included as special cases of the distributed loadings.

Boundary and Interface Conditions at the Stations
Since Equation (18) is a fourth order differential equation, four boundary or interface conditions for each segment have to be defined to yield a unique solution. According to Reference [39], the boundary conditions on the left end (station (1)) can be defined bỹ and on the right end (station (N)) bỹ whereŵ (1) prescribed harmonic moments at the left and right boundary, and f (•, ) is a linear function in terms of • and , which depends on the concentrated elements at the boundaries.
In case of classical boundary conditions (no concentrated elements at the boundaries), the displacementw (•) or shear forceQ (•) and rotationφ (•) or bending momentM (•) are prescribed. The most common types of classical boundary conditions are the clamped If concentrated elements are present at the boundaries, a coupling of the displacement and shear force and/or rotation and bending moment appears in the boundary conditions. The forces and moments due to concentrated elements acting on the beam boundaries are shown in Figure 3a,b (station (1) and station (N), respectively). Using Newton's second law of motion at the station (1) leads tõ and at station (N) results iñ whereF At an intermediate station (i), the displacement and rotation have to be continuous, which results inw where (i = ) and X − i and X + i are coordinates infinitesimal to the left and right of the station (i). The forces and moments at an intermediate station (i) are illustrated in Figure 4. According to Newton's second law of motion, the equilibrium of forces in z-direction and moments about the y-axis are given bỹ P point loadings, which are directly located at the station (i). The harmonic governing equation of an uniform Timoshenko beam given in Equation (18), combined with the boundary conditions in Equations (23)-(30) and the interface conditions in Equations (31)-(34), lead to a well-posed problem.

Numerical Assembly Technique
In the Numerical Assembly Technique (NAT), analytical solutions of each beam segment are used to fulfill the boundary and interface conditions. Since Equation (18) is a linear differential equation, the solution within a beam segment can be given bỹ as the homogeneous and particular solutions of the governing equation.

Homogeneous Solution of the Governing Equations
The general homogeneous solution of the harmonic governing equation (Equation (18)) is obtained by setting the external forces and moments to zero (q(x) = 0,m(x) = 0). Assuming a solution of the formw h (x ) = c w e j k x , and plugging it into Equation (18), results in the characteristic equation with the solutions and At the so-called cut-off angular frequency ω c , the characteristic of the third and fourth root changes; therefore, the form of the solutions is also different. For the undamped case, this frequency can be computed analytically, while, for the generally damped case, only a numerical approach is feasible. Using the roots in Equations (36) and (37), the final solutions for the four field-variables are given by with the constants In Equations (39) and (40), the field variables are gathered in the column vectorx h (x ) and the arbitrary constants c • in the column vectors c and c * . The state variable matrices B (x ) and B * (x ) below and above the cut-off frequency are slightly different, since the entries in the matrices are scaled to unity within the segment length (x = [0, L ]). This scaling has certain numerical advantages, since the terms with positive real parts in the exponents would grow very rapidly and lead to a bad conditioning of the resulting system matrices. In Section 3.3, the boundary and interface conditions are used to calculate the arbitrary constants c • .

Particular Solutions of the Governing Equations
The particular solutions of the governing equationw p (x ) fulfill the right-hand side of Equation (18). Depending on the type of loading, different solution procedures are necessary. In case of concentrated loadings, the Fourier transform [40] and the residue theorem and Jordan's Lemma [41] are applied to calculated the particular solutions. For generally distributed loadings, a two-step approach is required. First, the loading functions are approximated by a Fourier extension (continuation) method [33], and then the Green's function method [42] is used to compute the particular solutions. In the subsequent sections, an overview of the applied mathematical principle is given, and the particular solutions for point loadings and generally distributed loadings are derived.

Fourier Transform, Residue Theorem, and Jordan's Lemma
The Fourier transform is well suited to solve inhomogeneous differential equations, due to its operational property for derivatives [40] whereĝ(k) is the spatial Fourier transform F {•} of the function g(x) and n ∈ N 0 . The inverse Fourier transform F −1 {•} involves the evaluation of integrals on the real axis ranging from −∞ to ∞. A powerful tool to evaluate such integrals is the residue theorem [41]. The residue theorem leads to the integral formulas [41] ∞ −∞ĝ where z + i are poles in the upper half plane (total number s + ), z − i poles in the lower half plane (total number s − ), p i poles on the real axis (total number m), and Res{•} is the residue at the pole. The integral formulas in Equation (45) with | • | the absolute value [41]. Several methods are available to evaluate the residues. Since only simple poles appear in the following calculations, the formula [41] Res can be applied, which is valid for simple poles z i = ∞.

Concentrated Loads
Using a Fourier transform of Equation (18) with respect to the spatial coordinate x , the displacement of the beam segment in the transformed domain is given bŷ withŵ p (k),q(k), andm(k) as the spatial Fourier transforms ofw p (x),q(x), andm(x). The inverse Fourier transform leads to the particular solution in an integral form which can be analytically evaluated for concentrated loadings.

Point Force
A point force with the amplitudeF P acting at the coordinate x = X F is given bỹ with δ(•) the dirac delta function. The Fourier transform of the point force is given bŷ Plugging the Fourier transform of the point force into Equations (48) and (49) leads to the displacement and the evaluation of the integral with the residue theorem results iñ Due to the characteristic change of β at the cut-off frequency ω c , two different cases have to be considered. The upper signs in Equation (53) belong to ω ≤ ω c , and the lower signs to ω ≥ ω c . This compact representation of the two cases is adopted in all following equations. The remaining field variables are computed by Equations (19)- (21) and are given bỹ with sgn(•) the signum function. The particular solutions of the field variables are gathered in the column vectorx

Point Moment
The procedure for the derivation of the particular solutions due to a point moment is completely equal to the point force. For completeness, the solutions of the field variables are stated

Generally Distributed Loads
The particular solution functions for distributed loads are derived using the Green's function method, which uses the response due to a unit point source (Green's function) and an integration over the distributed load region [42]. In Figure 5a,b, the limits and coordinate systems used for distributed forces and moments are shown. Setting the point loadsF P andM P in Equations (53) and (57) to one, multiplying the solutions with the distributed loadingsq(x q ) andm(x m ), and integrating over the load region lead to the particular solutions which are given in a general integral form. Applying a change of variables to normalized local coordinates with respect to the limits of the distributed loadinḡ simplifies Equations (61) and (62) tõ with normalized angular wavenumbers.

Generally Distributed Force
Since the evaluation of the integral in Equation (65) strongly depends on the actual distribution of the excitation and has to be performed for any problem specific function individually, a different approach is taken. An approximation of the distributed force with the so-called Fourier extension (continuation) method is applied to allow for a systematic procedure to compute the particular solution for any given force distribution. In the Fourier extension method, a non-periodic function, which is defined on the interval [−1, 1], is approximated by a Fourier series that is periodic on an extended interval [−T, T] (T > 1) [43,44]. Therefore, the generally distributed force, given in normalized local coordinatesx q , is approximated by the Fourier series with T > 1 and n corresponding to the order of the approximation. A very efficient way to calculate the complex coefficients d qk is presented by Adcock and Huybrechs [43], which is called the numerical discrete Fourier extension with equispaced sampling points. Efficient and stable numerical algorithms of this method are developed in References [44,45]. A detailed description of the algorithm, which is applied in this paper, can be found in Reference [44], and additional information on the convergency rate and stability of the Fourier extension method is given in, e.g., References [43,46]. Plugging the Fourier extension of the distributed force (Equation (69)) into Equation (65), and analytically evaluating the integrals, leads to the particular solutioñ with the integrals I 1 (x , k) and I 3 (x , k) presented in Appendix A. The remaining field variables can either be computed by an integration of the point solutions (Equations (54)-(56)) or by applying Equations (19)- (21). The results for the rotation, bending moment, and shear force are given bỹ where the integrals I 2 (x , k) and I 4 (x , k) are also presented in Appendix A. The particular solutions for a generally distributed force are only quasi-analytical, since a small approximation error remains due to the Fourier extension method. This error can even be reduced to double-precision accuracy if the order n of the Fourier extension is high enough.

Generally Distributed Moment
The calculation of the particular solutions due to a generally distributed moment is completely equivalent. The distributed moment is approximated by a Fourier series and evaluating the resulting integrals in Equation (66) leads tõ

Assembly Procedure and Solution Process
The homogeneous and particular solutions, derived in Sections 3.1 and 3.2, are used to fulfill the boundary and interface conditions stated in Equations (23)- (34). This procedure leads to a system of linear equations of the form A c = b, with the 4M × 4M system matrix A and the 4M × 1 right-hand side vector b. The arbitrary constants c of each segment are gathered in the 4M × 1 vector c. Solving the system of linear equations for c results in a quasi-analytical solution of the beam vibration problem stated in Section 2.1. A detailed description of the assembly procedure can be found in the authors' previous paper [27].

Numerical and Experimental Validation Examples
In the following sections, numerical and experimental validation examples are presented. In Section 4.1, a general beam system is investigated using NAT and FEM, and the numerical solutions are compared with respect to accuracy and computational efficiency. Additionally, the behavior of different viscoelastic material models is examined. A validation of NAT and the fractional Zener material model with measurement data is shown in Section 4.2. Two different test setups with materials having low (steel) and high (PVC) internal damping are used to show the broad applicability of the fractional Zener model and the presented numerical technique.
All computations are performed on an Intel ® Xeon ® E3-1270 processor (4 × 3.6 GHz) with 32 GB RAM and a Windows 10 operating system. The software package MATLAB ® R2020b is used to implement NAT, while the FEM models are built in the commercial software Abaqus FEA ® 2017.

Numerical Validation Example
In Figure 6, a multiple-stepped beam with circular cross-sections including different loadings and support-conditions is shown. The beam with total length L = 1.5 m is divided by N = 7 stations into M = 6 uniform segments. Each segment of the beam has a constant circular cross-section, with the diameters d 1 = 0.08 m, d 2 = 0.1 m, d 3 = 0.12 m, d 4 = 0.15 m, d 5 = 0.11 m, and d 6 = 0.09 m. The material of the beam is uniform, and a polymer with high internal damping is chosen to show the dissipative behavior. In Weiß et al. [47], material parameters for PVC are presented, where the fractional Zener model is used to describe the viscoelastic material characteristics. The viscoelastic parameters, adapted to the notation in this paper, are a E 0 = 3.51 × 10 9 N m 2 , a E 1 = 8.5486 × 10 6 Ns α m 2 , b E 0 = 0.002106 s α , and α E = 0.6597 [47]. The material is isotropic with the Poisson's ratio ν = 0.38, and the density of the beam is ρ = 1450 kg m 3 , which leads to a total beam mass of 23.4 kg.  The beam is excited by a harmonic concentrated force with the amplitudeF P = 20 N at x = 0.2 m and a point moment with the amplitudeM P = 7 Nm at x = 0.53 m. Additionally, a distributed force acts in the fifth segment (X aq = 0 m and X bq = 0.15 m) given by the functionq(x 5 ) = 200 + 50 · sin(64 (1.05 + x 5 )) · (1.05 + x 5 ) 6 N m , resulting in a total force of 29.9 N.
The B22 element, which is a planer 3-node quadratic beam element having the Timoshenko theory implemented, is used to build the FEM models. To show the accuracy of the proposed computational technique, a FEM model with a very small element size of 0.001 m (6002 degrees of freedom) is applied, which guarantees a highly accurate reference solution. Furthermore, a coarse FEM model (70 degrees of freedom) with a similar computational time as NAT is used to illustrate the computational efficiency of the method. Since the commercial software package Abaqus FEA ® 2017 only supports a complex Young's modulus, the shear modulus is kept real-valued for the comparison with G = a E 0 2(1+ν) . For a detailed description of the finite element modeling process, the reader is referred to, e.g., Reference [48].
The frequency response with 2000 frequency steps in the interval of 0 Hz to 2000 Hz (10 modes) for the response point at x = 0.225 m is shown in Figure 7.  Figure 7, it is apparent that NAT gives highly accurate results, since the solutions of the reference model and NAT are not distinguishable. If no distributed loadings are present, NAT even leads to the analytical solution of the problem, due to the analytical homogeneous and particular solutions. The boundary and interface conditions are fulfilled with double-precision. For distributed loadings, a small error occurs when the actual loading is approximated by the Fourier extension method. In the given example, n = 32 terms are used in the approximation ofq(x 5 ), which leads to an excellent agreement. The coarse FEM model, which has a similar order of computational time as NAT, already suffers from certain approximation errors in the higher frequency range, as illustrated in Figure 7. This clearly shows the computational efficiency of NAT compared to element-based techniques. The total computational time of NAT for 2000 frequency points is 7.95 s, where 6.13 s are used to calculate the displacement, rotation, bending moment, and shear force at 1500 locations along the x-axis. Most of the computational time is, therefore, spend on the post-processing step, which can be reduced by a reduction of the evaluation points. Another advantage of NAT compared to FEM is the accuracy of derived quantities, such as the bending moment or shear force, since, in FEM, these quantities are approximated with a lower order polynomial.
In Figure 8, a comparison of the frequency-dependent loss factor η and storage modulus Re[E ] of the classical Zener model (standard linear solid) and the fractional Zener model is shown. The material parameters of Weiß et al. [47] are used, which are the result from a best fit to measurement data of a PVC beam and plate. While both models show a qualitatively similar behavior having one peak in the loss factor and a monotonically increasing storage modulus, the additional flexibility of the fractional model allows for a better control of the slope in the loss factor curve [11]. Therefore, the actual properties of the viscoelastic material (PVC) are captured more accurately. In Figure 9, the difference of both material models is also clearly visible in the frequency response for the beam system presented in Figure 6.  The difference in the storage modulus leads to a shifting of the eigenfrequencies between the models, while the different loss factors result in lower or higher amplitudes. The ninth mode at ≈1540 Hz and sixteenth mode at ≈3600 Hz are indicated by bullets and squares in Figures 8 and 9. Since the storage modulus and the loss factor of the classical Zener model are higher at the ninth mode, the eigenfrequency is slightly shifted to the right and the amplitude is reduced compared to the fractional model. At the sixteenth mode, the storage modulus of the classical model is higher, while the loss factor is lower compared to the fractional model. Therefore, the eigenfrequency and the amplitude are greater for the classical material model.

Comparison with Measurement Data
In this section, the presented numerical method and material model is compared to vibration measurements. Two different cases are considered. In the first case, a beam with a low material damping made of steel is investigated, while the second case shows a beam with high material damping (PVC).
The general test setup of the measurements is equal for both cases. An electrodynamic shaker (type: LDS V406, frequency range: 5 Hz to 9 kHz, maximum force: 196 N) is used to excite the system, and a force sensor (Brüel & Kjaer Type 8230-001, maximum force: 220 N, linearity error at full scale: ≤±1%) measures the excitation force. A triaxial acceleration sensor (Brüel & Kjaer Triaxial Type 4506, frequency range: 0.6 Hz to 3000 Hz) captures the response of the beam. The system is excited with a sinusoidal force at constant frequency, and the response is measured at steady-state conditions. The frequency step is 1 Hz, and both time signals are transformed into the frequency domain. The ratio of the beam displacement and force amplitude at the excitation frequency is calculated to get the Frequency Response Function (FRF) of the system.
In Figure 10a,b, the test setups of the steel and PVC beam are shown. The multiple-stepped beam a lies horizontally during the test, and fishing lines d are used to hold the beam in position. The force sensor including the shaker and stinger b and the acceleration sensor c are located in the same horizontal plane. Both sensors are mounted to the surface of the beam using a thin layer of wax. The fishing lines are used to get free boundary conditions at the ends of the beam. The NAT model used to compute the response of the test setup is illustrated in Figure 11. The beam has two steps along the x-axis with the diameters d 1 = d 4 = d 5 = 0.04 m and d 2 = d 3 = 0.05 m. The station parameters and locations are listed in Table 2. The excitation due to the electrodynamic shaker is modeled by a simple point force at station (3), since the introduction of a constant distributed force over the force sensor region has a negligible effect on the system response. The response point of the system is located at station (5). The masses of the sensors are included as lumped masses at the corresponding stations. The viscoelastic material parameters of the fractional Zener model are defined by a E 0 = 1.926 × 10 11 N m 2 , a E 1 = 6.887 × 10 9 Ns α m 2 , b E 0 = 0.035 s α , and α E = 0.3, which have been found through a simple optimization approach in MATLAB ® R2020b. The stated parameters are rather different compared to the values presented in Reference [35], which results from the low and nearly frequency independent damping of the material and the limited frequency range. Several parameter combinations lead to a very similar response in the viewed frequency range. Measurements up to a higher frequency would be required to resolve this issue, which is not feasible with the presented test setup. The material is isotropic with ν = 0.30, and the density is given by ρ = 7855 kg m 3 . The shear correction factor is set to k S = 0.8507, which is computed with the formulas for circular cross-sections shown in Reference [49].
In Figure 12, a comparison of the measured and calculated FRF is illustrated. There is a general good agreement of both results over the complete frequency range. At higher eigenfrequencies, the measurements are less reliable, since the limits of the applied measurement equipment is reach. Overall, the fractional Zener model is capable of representing the low, nearly frequency independent damping of steel without the drawbacks of the hysteric damping model. The parameters of the fractional Zener model used to describe the viscoelastic behavior of the PVC are given by a E 0 = 3.272 × 10 9 N m 2 , a E 1 = 8.565 × 10 6 Ns α m 2 , b E 0 = 0.00219 s α , and α E = 0.649, which are similar to the parameters presented in Reference [47]. The density of the material is given by ρ = 1377 kg m 3 , and an isotropic behavior with ν = 0.38 is assumed. The shear correction factor is calculated as in the previous case, which results in k S = 0.848.
A comparison of the measured and calculated FRF is illustrated in Figure 13. The results are in a good agreement over the complete frequency range up to 1000 Hz. The three eigenfrequencies at around 180 Hz, 430 Hz, and 870 Hz, as well as the amplitudes, are predicted very accurately by the NAT model. It is apparent that the internal damping of the PVC beam is considerable higher compared to the steel beam, since the peaks at the eigenfrequencies are less pronounced. The fractional Zener model is, therefore, also able to describe materials with high internal damping.

Conclusions
In this paper, a highly accurate and efficient numerical method for the dynamic analysis of viscoelastic beam systems, called Numerical Assembly Technique, has been presented. The viscoelastic behavior of the beam material is described by the fractional Zener model, which allows for the characterization of a wide class of materials with a low number of parameters. Analytical homogeneous solutions of the resulting governing equations have been presented, and analytical solutions for point force and moment excitation have been derived. Semi-analytical solutions for arbitrarily distributed loads have been calculated by an approximation with the Fourier extension method. A comparison to the Finite Element Method has shown the accuracy and efficiency of the proposed numerical method. Vibration measurements for a material with low damping (steel) and high damping (PVC) have been used to validate the material model, showing an overall good agreement between the calculated and measured results.

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:

Appendix A. Integrals
In this section, the evaluations of the integrals appearing in the particular solutions of generally distributed loadings are presented. Since the absolute value function and the signum function appear in the integrals, three different casesx ≥ 1,x ≤ −1 and |x | < 1 have to be considered during the evaluation. The final results are Tx − e −j( k π T +ᾱ (x +1)) j √ 2 T k π T +ᾱ + e j k π Tx − e j( k π T +ᾱ (x −1)) j √ 2 T k π T −ᾱ |x | < 1,