Steady-state linear harmonic vibrations of multiple-stepped Euler-Bernoulli beams under arbitrarily distributed loads carrying any number of concentrated elements

In this paper, general beam vibration problems with several attachments under arbitrarily distributed harmonic loading are solved. A multiple-stepped beam is modelled by the Euler-Bernoulli beam theory and an extension of an efficient numerical method called Numerical Assembly Technique (NAT) is used to calculate the steadystate harmonic response of the beam to an arbitrarily distributed force or moment loading. All classical boundary conditions are considered and several types of concentrated elements (springs, dampers, lumped masses and rotatory inertias) are included. Analytical solutions for point forces and moments and polynomially distributed loads are presented. The Fourier extension method is used to approximate generally distributed loads, which is very efficient for non-periodic loadings, since the method is not suffering from the Gibbs phenomenon compared to a Fourier series expansion. The Numerical Assembly Technique is extended to include distributed external loadings and a modified formulation of the solution functions is used to enhance the stability of the method at higher frequencies. The method can take distributed loads into account without the need for a modal expansion of the load, which increases the computational efficiency. A numerical example shows the efficiency and accuracy of the proposed method in comparison to the Finite Element Method. c © 2020 University of West Bohemia. All rights reserved.


Introduction
In engineering applications, real structures can be frequently modelled by an assembly of beams with attached concentrated elements like lumped masses, springs and dampers, e.g. lever arms, shafts, aeronautical structures, robotic arms etc. The vibration characteristics of a beam system are of major interest to engineers, since even small dynamic loadings can lead to catastrophic failure of machines or structures [28].
Several analytical and semi-analytical methods are available to analyze the free vibrations of beams. The Numerical Assembly Technique (NAT) [9][10][11]18,19,[31][32][33][34][35]42,43,47,[49][50][51]53] uses the analytical solutions of the homogeneous governing equation within one uniform beam segment and enforces the boundary and interface conditions. In this method, the complete beam has to be divided into several segments if concentrated elements or changes in the material parameters or cross-section appear. Therefore, the size of the system matrix rises with the complexity of the beam vibration problem, but the implementation is straight forward and the system matrix has a banded structure.
In the Transfer Matrix Method (TMM) [5,6,[44][45][46], the beam state (displacement, rotation, bending moment and shear force) on one end of a uniform beam segment is transferred to the other end by a transfer matrix. Therefore, the system matrix for two-dimensional beam vibration problems has a constant size of 2 × 2 independent of the number of concentrated elements or changes in the material parameters or cross-section, since unknowns in the intermediate beam segments are eliminated.
The Dynamic Stiffness Method (DSM) [4,23] applies a dynamic stiffness matrix, which connects the forces and moments within a beam segment with its displacements and rotations. The complete system matrix is built by an assembly of the exact dynamic stiffness matrices of each beam segment, similar to the assembly process in the Finite Element Method (FEM). Therefore, the method is very general and leads to analytical results for the free vibration problem of beams.
In the Green function method (GFM) [24][25][26][27]41], the analytical responses of an uniform beam segment to concentrated loads (Green functions) are used to include concentrated elements in the model. The remaining unknowns are the displacements and rotations of the beam at the positions of the concentrated elements. Therefore, the size of the system matrix also increases with the complexity of the beam vibration problem, but slower compared to NAT.
A very similar method is the Generalized Function Approach (GFA) [8,15,17], which uses generalized functions (Dirac delta functions, Heaviside step functions) to treat the discontinuities in the response of the beam due to concentrated elements and changes in the material parameters and cross-section. Compared to GFM, the GFA leads to a 2 × 2 system matrix for the twodimensional beam vibration problem independent of the number of concentrated elements and changes in the material parameters or cross-section. Therefore, the method is computationally more efficient.
All these methods lead to a frequency-dependent system matrix and the eigenfrequencies and associated mode shapes are found by setting the determinant of the system matrix to zero. This results in a transcendental equation, which has to be solved numerically.
Generally, all mentioned methods can be used to analyze forced vibrations of beams by applying the modal superposition method (MSM) [40]. There, the response of the beam to an arbitrary external load is built by an infinite series of the orthogonal mode shapes, which is truncated at a certain term. Due to the truncation, only approximated results are obtained and generally, the procedure is computationally demanding. GFA and MSM have been applied in [15] to find the forced response of an Euler-Bernoulli beam with an arbitrary number of Kelvin-Voigt viscoelastic attachments and in [17] for plane beam structures with mass-spring subsystems and rotational joints.
If the structure is excited by point forces or moments, the external loadings can be included in the interface conditions between two beam segments. This approach has been used in [36] to calculate the forced response of a beam with multiple supports applying TMM. Similar, the response of a beam with several concentrated attachments excited by a harmonic point force is determined by NAT in [30,48,52]. An additional splitting of the beam is necessary in this approach to include the external loadings, which leads to a larger system matrix and therefore, a higher computational time. Furthermore, only concentrated loads can be treated and an extension to distributed loads is not feasible.
A very common approach to find the forced vibration response of a beam is the calculation of the so-called Green's function. The Green's function is the response of the beam to an arbitrarily positioned unit point load. Due to the linearity of the beam vibration problem, the response to a distributed load can be found by integrating the Green's function multiplied by the load distribution [22]. This approach has been presented in [1,12,20,29] (GFM) and [7,16] (GFA) for different beam vibration problems with and without attachments. One of the main challenges of this approach is the integration of the Green's function, which can be time consuming and is not always possible in an analytical way.
In this paper, the Numerical Assembly Technique is extended to allow for the analysis of general beam vibration problems with any number of concentrated attachments (translational and rotational springs, translational and rotational dampers, lumped masses and rotatory inertia) excited by an arbitrarily distributed harmonic loading. Therefore, the analytical solutions of infinitely extended Euler-Bernoulli beams excited by point forces and moments and polynomially distributed forces and moments are derived in closed-form. If the external loading is arbitrarily distributed, an approximation of the load by the Fourier extension method is applied, which leads to a representation of the load as partial sum of complex exponential functions. The response of an infinitely extended Euler-Bernoulli beam to such a partial sum is also analytically calculated in closed-form. While a classical Fourier series expansion would suffer from the Gibbs phenomenon if the approximated function is non-periodic, the Fourier extension method [2,3,37] overcomes this limitation. Therefore, the convergency rate of the Fourier extension method is, in general, very high and an efficient implementation of the method is also available [37].
Compared to the standard NAT, a different solution of the homogeneous governing equation is applied, which leads to a better conditioning of the system matrix, especially for higher frequencies and long and thin beams. This enhances the stability and accuracy of the method. Furthermore, an additional splitting of the beam into sub-segments, if an external load appears, is not necessary in the proposed extension of NAT, which leads to a smaller system matrix and higher computational efficiency.
The proposed method is quasi-analytical in the sense that the governing equations are fulfilled exactly and only minor errors are introduced in the enforcement of the boundary and interface conditions due to double-precision floating point arithmetic when solving the resulting system of linear equations. A small approximation error also occurs if the external loading is not concentrated or polynomially distributed, but the results can be made arbitrarily accurate by increasing the order of the Fourier extension when approximating the generally distributed load.

Problem description and beam theory
In this section, a general two-dimensional beam vibration problem with several concentrated elements and external loading is described and the Euler-Bernoulli theory for the transverse vibration of straight and uniform beams is outlined.

General beam vibration problem with concentrated elements and external loading
In Fig. 1, a general beam vibration problem with concentrated elements and external loading is shown. The beam with total length L is divided by (N) stations into M = (N − 1) segments. The first (1) and last station (N) are located at the beam boundaries (x = 0 and x = L) and additional stations (i) have to be included if concentrated elements (translational springs k r , lumped masses m (i) or rotatory inertia Θ (i) ) are present or a change in the beam material or geometric properties appears. The location Each beam segment is uniform with the length L = (X i+1 − X i ) (i = ), constant Young's modulus E , density ρ , cross-section area A and second moment of area about the y-axis I . For brevity, it is assumed that all types of concentrated elements and a change in material and geometric properties appear simultaneously at each station. Simply setting certain parameters of the concentrated elements to zero or keeping the material or geometric properties unchanged, allows for every possible configuration at the stations.
The beam is loaded by point forces F P (t) at the positions X F , point moments M P (t) at the positions X M , distributed forces q(x, t) and distributed moments m(x, t). Since steadystate harmonic vibrations of the beam are investigated, the external loads are assumed to vary harmonically with angular frequency ω and, therefore, are given by (1) with t the time, j = √ −1 the imaginary unit andq(x),m(x),F P andM P the complex amplitudes of the external loads. If external loads with different angular frequencies are acting on the beam, the superposition principle can be applied, since the beam vibration problem is linear. For periodic excitation, the Fourier series representation can be used to transform the loading into its harmonic components.

Euler-Bernoulli beam theory
Each segment of the beam is modeled by the well known Euler-Bernoulli beam theory with the harmonic governing equation [21] where withw (x) the complex amplitude of the displacement in z-direction and δ(•) the Dirac delta function. In Eq. (2), harmonic vibrations at angular frequency ω are assumed and the timedependent term e jωt is dropped.

Boundary and interface conditions at the stations
The harmonic governing equation of the beam segment , given in Eq. (2), requires certain boundary and interface conditions to yield a unique solution. Apart from the displacement amplitudew (x), the field variables withφ (x),M (x) andQ (x) the complex amplitudes of the beam rotation, the bending moment about the y-axis and the shear force in z-direction, are introduced in the boundary and interface conditions. According to [40], the classical boundary conditions on the left end (station (1)) can be defined byw and on the right end (station (N)) bȳ whereŵ (1) andŵ (N ) are prescribed displacements,φ (1) andφ (N ) prescribed rotations,F If concentrated elements are present at the boundaries, the displacement and shear force and/or rotation and bending moment are coupled, which leads tō at the station (1) and results in at station (N). At an intermediate station (i) the continuities of displacement and rotation with (i = ) have to be fulfilled. In Eqs. (15) and (16) the locations X − i and X + i are infinitesimal to the left and right of the station (i). The equilibrium of forces in z-direction and moments about the y-axis at an intermediate station (i) results in P point forces and moments, which are located at the station (i). The harmonic governing equation of an uniform Euler-Bernoulli beam given in Eq. (2), the boundary conditions in Eqs. (7)- (14) and the interface conditions in Eqs. (15)-(18) form a well-posed problem, which can be uniquely solved by analytical or numerical methods.

Homogeneous solution of the harmonic governing equation
The general homogeneous solution of the harmonic governing equation (2) is obtained by setting the external loads to zero (q(x) = 0,m(x) = 0,F P = 0 andM P = 0). Assuming a solution of the formw h (x ) = c e j k x leads to the characteristic equation with the four roots k 1 = λ , k 2 = −λ , k 3 = −j λ and k 4 = j λ . Therefore, the general solution of the homogeneous governing equation in local coordinates x is given bȳ with c 1 , c 2 , c 3 and c 4 arbitrary constants. The lower index • h indicates the homogeneous solution of the differential equation. Generally, the solution in Eq. (20) is rewritten in terms of sin(•), cos(•), sinh(•) and cosh(•). In this paper, a slightly different approach is used by rewriting Eq. (20) as withc 1 ,c 2 ,c 3 andc 4 different arbitrary constants. Using the solution given in Eq. (21) guaranties that the amplitude of each function term remains smaller or equal to 1 within the segment span (0 ≤ x ≤ L ), which is advantageous for the numerical stability of NAT. A complete description of the state within a beam segment is given, if not only the displacement, but also the rotation, bending moment and shear force are known. If the state variables and arbitrary constants are gathered in the column vectors the state within a beam segment can be defined by a compact matrix equation as with the state variable matrix where the upper index • T denotes the transpose of a vector or matrix and B = E I , C = cos(λ x ), S = sin(λ x ), E + = e λ (x −L ) , E − = e −λ x . The relations in Eqs. (4)- (6) have been used to derive Eqs. (22) and (23).

Particular solutions of the harmonic governing equation
The total displacementw (x ) of a beam segment is given by the sum of the homogenous solutionw h (x ), derived in Section 3, and a particular solutionw p (x ), which fulfils the right hand side of the harmonic governing equation (2). In this section, the particular solution functions for different types of excitation are derived, using the Fourier transform [14], the residue theorem and Jordan's lemma [38] and the Green's function method [22].

Point forces and moments
The particular solution functionw p (x) for a beam segment under point force excitation can be derived by applying the Fourier transform to Eq. (2), which finally leads to the integral which can be evaluated by the residue theorem stated in [38]. This leads to the displacement due to an external point force excitation and the rotation, bending moment and shear force, computed through Eqs. (4)- (6), are gathered in a column vectorx The results for a point moment excitation can be derived in a similar way and are given bȳ with sgn (•) the signum function.

Polynomially distributed forces and moments
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 [22]. In Fig. 2, the limits and coordinate system used for distributed forces are shown. In this section, the distributed loadsq(x) andm(x) are represented by polynomials, which are defined in local coordinates (x q and x m ) by c qk x k q with X Aq ≤ x q ≤ X Bq (27) andm where N q and N m are the order of the polynomials, 0 ≤ X Aq < X Bq ≤ L and 0 ≤ X Am < X Bm ≤ L the limits of the polynomially distributed force and moment and c qk and c mk the coefficients of the polynomials. Using Eq. (25) withF p = 1 and the Green's function method leads to the particular solution for a distributed force, which is given by the integral The evaluation of the integral appearing in Eq. (29) is shown in Appendix A. The entries of the state variable vectorx p (x ) = w p (x ),φ p (x ),M p (x ),Q p (x ) T for the polynomially distributed force can be defined by Eq. (53) and for the polynomially distributed moment by Eq. (55). Therefore, the particular solution functions for polynomially distributed loads are given in an exact, analytical and closed form.

Generally distributed forces and moments
In general, the distributed loadsq(x q ) orm(x m ) may not be represented by polynomials, but by any arbitrary function. Plugging an arbitraryq(x q ) into Eqs. (29) leads to the particular solutions in an integral form. The evaluation of the resulting integrals can be time consuming, because an analytical integration is not always feasible and numerical integration schemes have to be applied. In this paper, a different approach is used, by approximating the generally distributed load by the Fourier extension (continuation) method [2,3,37]. The approximation of periodic functions is efficiently and accurately possible by a classical Fourier series, but the so-called Gibbs phenomenon prevents a fast convergence for non-periodic functions [37]. The Fourier extension method overcomes this shortcoming, while retaining most of the advantages of the classical Fourier series, which allows for an analytical closed-form particular solution for a generally distributed loading.
The limits and coordinate systems used for the generally distributed loads are identical to the polynomially distributed ones and are shown in Fig. 2. Using the transformations of variables to normalized local coordinates simplifies the integral in Eq. (29) tō In Eq. (32) the normalized angular wavenumberλ is given bȳ In the Fourier extension (continuation) 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) [3,37]. Therefore, the generally distributed force and moment, given in normalized local coordinatesx q orx m , are approximated by the Fourier series Tx q (34) andm Tx m (35) with T > 1 and n corresponding to the order of the approximation. Several different approaches exist to calculate the complex coefficients d qk and d mk , but the numerical discrete Fourier extension with equispaced sampling points [3] seems to be the most suited one for the problem at hand, since efficient and stable numerical algorithms are available [37]. The reader is referred to [37] for a detailed description of the numerical algorithm used in this paper and to [2] and [3] for an in-depth analysis of the convergence rate and stability of the

Numerical Assembly Technique
In the Numerical Assembly Technique (NAT), the homogenous solutionx h (x ), derived in Section 3, and the particular solutionsx p (x ), given in Section 4, are used to fulfil the boundary and interface conditions, stated in Section 2.3, and therefore, the harmonic governing equation (2) is exactly satisfied. Plugging the total solutionsx (x ) =x h (x ) +x p (x ) of each beam segment into the boundary and interface conditions leads to a system of linear equations Ac = b with c = [c 1 , · · · ,c , · · · ,c M ] T the arbitrary constants of the homogenous solutions. The system matrix A depends only on the type of boundary conditions and concentrated elements at the stations, while the right-hand side vector b is additionally affected by the prescribed boundary values and external loading.
If no concentrated elements are present at the boundaries, the four different classical types of boundary conditions (7)-(10) can be written in matrix form using the total solutionsx 1 (0) andx M (L M ): • Prescribed displacement in z-direction and rotation about the y-axis: • Prescribed bending moment about the y-axis and shear force in z-direction: • Prescribed displacement in z-direction and bending moment about the y-axis: • Prescribed rotation about the y-axis and shear force in z-direction: In Eqs. (37)-(40), the 2 × 4 matrices A (1) and A (N ) and the 2 × 1 right-hand side vectors b (1) and b (N ) are associated with the left (station (1)) and right boundary (station (N)). The left upper index in • B indicates the • th row in the state variable matrix B .
If concentrated elements are present at station (1) and (N), Eqs. (11) and (12) and Eqs. (13) and (14) have to be satisfied. Using the total solutionx 1 (0) of the first beam segment in Eqs. (11) and (12) leads to for the left station (1) and applying the solutionx M (L M ) of the last beam segment in Eqs. (13) and (14) for the right station (N). In Eqs. (41) and (42), the parameters of the translational concentrated elements at station (•) are gathered in g The 4M × 4M system matrix A and the 4M × 1 right-hand side vector b are an assembly of all sub-matrices A (i) and sub-vectors b (i) , where i = (1, . . . , N). The solution of the system of linear equations Ac = b leads to the unknown constantsc and therefore, the displacement in z-direction, the rotation about the y-axis, the bending moment about the y-axis and the shear force in z-direction are uniquely defined for the whole beam.

Numerical results and discussion
In this section, a beam vibration problem with general support-conditions and loadings is analyzed. The computational efficiency and accuracy of NAT for steady-state harmonic beam vibration problems are demonstrated. NAT is implemented in MATLAB R R2019a and the reference solutions are calculated with the commercial FEM software package Abaqus FEA R 2017. The planar and cubic beam element B23, which has the Euler-Bernoulli beam theory implemented, is applied in all FEM calculations and a mesh size of 0.001 m is used to guarantee accurate reference solutions. The system of linear equations in NAT are solved by the Matlab built-in function linsolve, which uses a LU-factorization with partial pivoting. All computations are performed on an Intel R Core TM i7 system (6 × 3.7 GHz, 32 GB RAM), which uses a Windows 10 operating system.  Table 1. The total beam length is 1.25 m and the total mass of the beam is 12.89 kg.
The four-step beam is shown in Fig. 3, with the locations of the six stations, the concentrated forceF P = 500 N at x = 0.2 m, the concentrated momentM P = 100 N m at x = 0.525 m, the polynomially distributed forceq 1 (x 2 ) = 400 (100 x 2 2 − 20 x 2 + 4) (total load 266.66 N) and the generally distributed forceq 2 (x 4 ) = 2 000 exp (sin 2 (50 x 4 )) (total load ≈ 687 N). Since springs and dampers are present at the first and last station, no classical boundary conditions have to be prescribed. The generally distributed forceq 2 (x 4 ) is approximated by the Fourier extension method with n = 32 for the calculations with NAT. Even though, the amplitudes of the excitations are all purely real-valued, the response amplitudes are complex-valued due to the damping in the system, which introduces a phase shift. The real part of the displacement amplitudes for an excitation frequency f = 2 100 Hz calculated with NAT and FEM is shown in Fig. 4a and the imaginary part in Fig. 5a. It is apparent that the results of NAT and FEM are practically identical, which can also be seen by the relative errors plotted in Figs. 4b and 5b.
The frequency response function for the four-step beam at x = 1.00 m from 10 Hz to 5 000 Hz, computed with NAT and FEM at 501 steps, is illustrated in Fig. 6. The amplitude (Fig. 6a) and phase (Fig. 6b) calculated with NAT and FEM are perfectly coincident over the whole frequency range, which shows the accuracy of NAT for the most general steady-state harmonic beam vibration problem stated in Section 2.1.
Only the first eigenfrequency can be clearly identified by the peak in Fig. 6a, since the higher modes are strongly damped. The general phase shift of the displacement amplitude to the excitation varies between −180 • and 180 • due to the local dampers in the stations. The total time to compute the frequency response function at 501 frequency steps with NAT and evaluate the displacement, rotation, bending moment and shear force at 125 equally spaced response points is less than 1.3 s, where 0.8 s are required for the post-processing steps. The computation efficiency of NAT is hardly effected by the additional damping, even though the system matrix becomes complex-valued.

Conclusion
The numerical example shows the overall excellent agreement of the results calculated with NAT and FEM. The results of NAT are expected to be highly accurate, since NAT is quasi-analytical in the sense that analytical homogeneous solutions and analytical particular solutions for point and polynomially distributed excitations are used. The only errors introduced in the solution process are the approximation errors of the generally distributed loads by the Fourier extension method and the small numerical error due to the double-precision floating point arithmetic when solving the system of linear equations.
The total computational time of NAT in the example is very low, which shows the high efficiency of the method for steady-state harmonic beam vibration problems. The computation time is mainly effected by the type of excitation, since the evaluation of the upper incomplete gamma function (polynomially distributed load) is rather time consuming compared to trigonometric functions (all other load types). Additionally, the total number of segments has a major impact on the computation time, since the size of the system matrix becomes larger with an increasing number of segments. A main advantage of the proposed extension of NAT is the fact that new stations have to be only included if concentrated elements or changes in the cross-section or material parameters appear. This is in contrary to similar methods, which require additional stations when external loading is present.
The system matrix in NAT is purely real-valued as long as no damping is introduced into the system and becomes complex-valued otherwise. Although, the system matrix becomes complexvalued, the computational efficiency of NAT is hardly effected by the presents of damping. It is also possible to add material damping, like a Kelvin-Voigt or Maxwell material model [13], by simply replacing the Young's modulus E by a complex dynamic modulus E * .
It needs to be noted that distributed moments is not included in the numerical results due to a lack of reference solutions, since no implementation of distributed moments are available in the commercial FEM software packages. There is no difference in the solution process between distributed forces and moments and therefore, accurate results are also expected for distributed moment excitation.
In general, the proposed method can be extended to any linear structural vibration problem as long as the governing equations are reducible to a system of ordinary differential equations. Therefore, NAT is also applicable for, e.g., the more accurate Rayleigh beam theory or the Timoshenko beam theory and coupled three-dimensional beam vibration problems.