On the numerical properties of high‐order spectral (Euler‐Bernoulli) beam elements

In this paper, the numerical properties of a recently developed high‐order Spectral Euler‐Bernoulli Beam Element (SBE) featuring a C1‐continuous approximation of the displacement field are assessed. The C1‐continuous shape functions are based on two main ingredients, which are an Hermitian interpolation scheme and the use of Gauß‐Lobatto‐Legendre (GLL) points. Employing GLL‐points does not only avoid Runge oscillations, but also yields a diagonal mass matrix when exploiting the nodal quadrature technique as a mass lumping scheme. Especially in high‐frequency transient analyses, where often explicit time integration schemes are utilized, having a diagonal mass matrix is an attractive property of the proposed element formulation. This is, however, achieved at the cost of an under‐integration of the mass matrix. Therefore, a special focus of this paper is placed on the evaluation of the numerical properties, such as the conditioning of the element matrices and the attainable rates of convergence (ROCs). To this end, the numerical behavior of the SBEs is comprehensively analyzed by means of selected benchmark examples. In a nutshell, the obtained results demonstrate that the element yields good accuracy in combination with an increased efficiency for structural dynamics exploiting the diagonal structure of the mass matrix.


INTRODUCTION
Despite the ever increasing computational power, beam and shell elements are still a vital part of structural engineering analysis and will remain so for the foreseeable future. Therefore, a continued research effort is directed at developing structural finite elements with certain properties tailored to specific applications. In the element libraries of commercial finite element (FE) tools, structural finite elements, including beam [1], plate [2,3], and shell [4,5] elements, play an integral part to allow for an efficient modeling of frame structures and thin-walled designs. Different from solid elements, many structural finite elements consider rotational degrees of freedom at the nodes, which enables them to efficiently model structures, where bending deformations are dominant. In the static analysis of frame structures, for example, it is usually preferable, but not always possible, that only one element is used per structural member [6], which means that special structural finite elements must be exploited. However, for high-frequency dynamics it is inevitable to employ several structural elements per member as the wavelengths of the propagating waves need to be accounted for. This is still numerically much more efficient compared to a discretization based on solid elements, where both a minimum amount of elements along the length and over the thickness of the structural member are required to accurately represent the bending part of the deformation field. This will, however, drastically increase the computational costs. As a result, despite of the rapid development of modern computers, nowadays structural elements are still recommended in the finite element analysis (FEA) of many structures, including bridges, piles, tunnels, fibers, and others. Therefore, a vast abundance of research directed at FE modeling with structural finite elements can be found in the wide body of literature [7][8][9][10]. In this work, we mainly focus on beam elements as a representative structural finite element, although the idea can be readily extended to plate and shell elements for more complicated problems, especially if cut Bogner-Fox-Schmit plate/shell elements, as introduced in Refs. [11,12], are employed.
Transient wave propagation in solids and structures is encountered in several applications such as non-destructive evaluation (NDE) [13,14], structural health monitoring (SHM) [15][16][17][18], impact response [19], seismology [20], acoustics [21][22][23], and biomechanics [24] to mention just a few. In this context, ultrasonic guided waves provide a robust tool for damage detection in thin-walled flat and curved structures, such as aircraft fuselages, ship hulls, pressure vessels, and pipelines [15,16,[25][26][27]. This is due to the fact that they can detect both surface and embedded damages of small size and can do so over a long distance at a very short time. Hence, they play an important role in both offline NDE and the emerging online SHM technologies. Note that numerical simulations of guide waves are often required to effectively design SHM systems because of the high costs of experimental investigations. Further, many SHM strategies employ model-based techniques or physics-informed data-driven models, which require the solution of inverse problems. In all these cases, to obtain the final solution many wave simulations must be performed, which results in high computational costs and hence, computationally efficient algorithms as developed in this contribution are of utmost importance.
Moreover, the waves for probing a structure need to have small wavelengths to detect small size damage and therefore, they must be excited at high frequencies (typically 50 to 1000 kHz). The classical finite element method (FEM) is usually adopted for performing such dynamic analyses. However, for high-frequency wave propagation problems, it takes a large numerical effort (or high computational time) because of the requirements of fine spatial as well as temporal discretizations [28,29] and it is prone to numerical dispersion and dissipation errors [30]. Thus, the classical FEM may be too inefficient for SHM applications. This has resulted in numerous efforts to develop improved methods for the accurate and efficient simulation of wave propagation in structures; comprehensive reviews of these activities can be found in Refs. [31,32].
The semi-analytical finite element (SAFE) method [33][34][35][36] and the wave finite element (WFE) method [37,38] are representatives of methods which are efficient, but restricted by the assumptions of either uniform or periodic cross-sections. Furthermore, these methods are computed in the frequency domain and thus, require an inverse Fourier transform to recover the time history response. This process may be prone to wrap-around errors, costly, and even unsteady in practice [39].
On the other hand, the frequency domain spectral element method (FDSEM) proposed by Doyle [39] eliminates the restriction of uniform or periodic cross-sections. In this method, the governing partial differential equations are transformed into the frequency domain using either Fourier, Laplace, or wavelet transforms, and then solved analytically to obtain a frequency-dependent dynamic stiffness matrix [40]. The method has been employed for wave propagation analysis of isotropic and composite straight, or curved beams [41][42][43][44][45] and plates having an infinite length or Levy-type boundary conditions (BCs) along one direction [46,47]. Moreover, guided wave propagation in laminated and delaminated beams has been reported in Refs. [48,49]. Aside from the numerical challenges associated with the inverse transforms, particularly when the frequency sampling is dense, the FDSEM is challenging to use in two-dimensional waveguides of complex geometries and BCs.
Another approach, called enriched FEM, has been recently proposed for solving general wave propagation problems, in which the polynomial interpolation functions of the classical FEM are enriched with local element-domain wave packet functions satisfying the partition of unity condition [50][51][52]. It has been shown to be accurate and computationally efficient for a large variety of wave propagation problems in one-and two-dimensional (planar) as well as axisymmetric domains [50][51][52]. However, the method is yet to be extended to beam and plate/shell models.
The time-domain spectral element method (SEM), proposed by Patera [53], is possibly the most widely used technique for solving guided wave propagation problems. The reason being that it directly provides the time signals (without requiring any transformation) and fully retains the advantages of the classical FEM for modeling complex geometries. It employs high-order polynomial interpolation functions, but instead of equispaced nodes, as used in the classical FEM, the nodes are located at the roots of orthogonal polynomials [54,55]. This change helps to alleviate large oscillations in the shape functions near element boundaries, known as Runge's phenomenon [56], and also yields spectral convergence. Additionally, the nodal quadrature technique may be utilized to diagonalize the mass matrix, increasing the efficiency of explicit time integration schemes. The error associated with the use of GLL points has been analyzed in Ref. [57]. Certain numerical issues of wave propagation modeled by the SEM are discussed in Ref. [58], in which the difference between equidistant points, GLL points, and Gauß-Lobatto-Chebyshev (GLC) points are discussed. In Ref. [54], a rod is modeled using high-order spectral elements (SEs) with an additional concentrated mass acting as damage for wave propagation analysis. It is concluded that the location of the mass can be successfully detected when comparing with experimental data. In Ref. [59], an accurate layerwise zigzag theory has been developed for the bending analysis of laminated beams. In a nutshell, it can be stated that SEs based on the theory of elasticity have been derived for the analysis of Lamb wave propagation in 1D [54], 2D, and 3D [60][61][62][63] domains and thus, are readily available for different areas of application. Due to the fact that all elements discussed above are based on the theory of elasticity, C 0 -continuous shape functions are sufficient to fulfill the requirements of the weak form of the underlying governing equations.
For thin-walled structures, it is important to use beam and shell theories that take advantage of their small thickness and make assumptions on the variations of the displacement field along the thickness direction to gain a computational advantage. Accordingly, SEs have been developed for singly-layer and laminated beams and plates based on the firstorder shear deformation theory (FSDT) [55,[64][65][66][67] and higher-order layerwise theories [68][69][70]. The application of SEM in detecting breathing cracks using guided waves is, for example, reported in Refs. [71,72]. These elements, too, are based on C 0 -continuous interpolations. Note that there are several beam and shell theories such as the classical beam/shell theories, the advanced third-order shear deformation theory (TOT) [2], and the efficient layerwise zigzag theory [43] that require C 1 -continuity for the deflection interpolation. Although, the refined third-order beam theory developed by Reddy [6] features a more realistic shear stress distribution along the cross-section of the beam and zigzag theories, such as the one developed by Cho and Parmerter [73], yield very accurate results, until recently, no SEs based on such theories have been available. This gap was filled by Kapuria and Jain [74] by presenting a SE for isotropic Euler-Bernoulli beams introducing C 1 -continuous Hermite-type interpolation functions. Using this interpolation technique, they have further developed SEs based on the refined TOT for isotropic [75] and laminated composite beams [76]. Composite materials, especially fiberreinforced ones, have drawn considerable attention in recent years due to their excellent mechanical properties [77,78] and are thus, highly researched. In fiber-reinforced composites (FRCs), fibers are sometimes modeled with bar elements [79], which fail to account for the bending stiffness of the fibers and hence, the application of beams elements is often required [80,81], which marks a future area of application of the beam elements investigated in this article.
The high-order C 1 -continuous spectral Euler-Bernoulli beam elements (SBEs) proposed by Kapuria and Jain in Ref. [74] exhibit several properties, which make them ideal for the analysis of transient problems. It has been shown that compared to low-order formulations the developed SBEs yields very accurate results and a significantly improved convergence for the free vibration response of beams, which is expressed by drastic savings in terms of degrees of freedom or computational time. Moreover, all theories of thin-walled structures requiring a C 1 -continuous representation of the displacement field can be readily implemented in conjunction with the developed SBEs, making them ideal for different types of problems as demonstrated in Refs. [74][75][76]. Due to the excellent properties of the SBEs, a detailed analysis of the numerical properties of this class of elements is of utmost importance, but unfortunately still missing. Therefore, the overarching goal of this paper is to fill this knowledge gap. This analysis paves the way for a wide-spread application of this element type to various problem classes.
As discussed in the previous paragraph, the main novelty lies in a comprehensive analysis of the numerical performance of the SBEs, which has not been conducted yet. To this end, we recall that the proposed SBEs are based on two main ingredients: (i) the use of an Hermitian interpolation scheme and (ii) the application of GLL points, which does not only avoid Runge oscillations, but also yields a diagonal mass matrix when exploiting the nodal quadrature technique. It is well-known that a diagonal mass matrix is a prerequisite for the efficient analysis of problems in high-frequency dynamics and therefore, it is of special interest to investigate the ramifications of using a lumped mass matrix, which is obtained at the cost of an under-integration of the inertia term in the weak form. Accordingly, a special focus of this paper is placed on the assessment of the effect of the (numerical) integration error on properties of the SBEs, such as the condition number of the element matrices, the attainable rates of convergence (ROCs), and so forth. To this end, the numerical behavior of the SBEs is comprehensively analyzed through selected benchmark examples, revealing pertinent numerical properties of this element type.

SPECTRAL BEAM ELEMENTS
SEs commonly rely on a tensor-product formulation of Lagrangian interpolation polynomials defined on non-equidistant nodal grids [53,[82][83][84][85]. One feature, they all have in-common, is that the derived shape functions are only C 0 -continuous 1 .
This, however, does not suffice when dealing with beam and shell elements based on the classical lamination theory (CLT) [46], a refined TOT [2,6], or layerwise zigzag theories [43,59]. These approaches require a C 1 -continuous displacement field, which can be readily achieved by either exploiting the isogeometric analysis (IGA) paradigm [87] or by considering an Hermitian interpolation technique. The latter approach is favored due to its simplicity and similarity to classical finite element procedures. In previous works, this idea has already been successfully implemented by Kapuria and Jain [74] for the purpose of wave propagation analysis in the context of SHM applications. For the sake of a self-contained description of the theory, the main ingredients needed to derive these C 1 -continuous SBEs are recalled in the remainder of this section.

Hermitian shape functions
Consider a simple Euler-Bernoulli beam element of length in its local reference frame (see Figure 1), featuring an arbitrary number of nodes N at the positions , where each node is associated with two degrees of freedom (DOFs), that is, and (= d ∕d ), corresponding to the transverse deflection and slope (rotation) at this node. In the current analysis, we assume that the beam axis is straight and therefore, the local element coordinate is mapped to the global coordinate using a linear relation Thus, the derivative of an entity with respect to the global coordinate □ , can be straightforwardly expressed by its derivative with respect to the local coordinate □ , and the length of the beam element Our main goal is to achieve a C 1 -continuous interpolation polynomial, which necessitates that both the deflection and slope must be prescribed at each node. Consequently, a beam element featuring N nodes requires shape functions of order = 2 N − 1. In a generic way, this interpolation function is written as Additionally, we also require the first (global) derivative of the displacement field, which yields the variation of the rotation along the beam axis By introducing the vector of all monomials of the ansatz and the vector of unknowns , we can express Equations (3) and (4) in matrix form as F I G U R E 1 -Node spectral beam element in its local coordinate system. with = [1 2 ⋯ 2 N −1 ], , = 2∕ [0 1 2 3 2 ⋯ (2 N − 1) 2 N −2 ], and = [ 0 1 2 ⋯ 2 N −1 ] T . Note that the unknowns of the ansatz do not hold a physical meaning. Therefore, it is preferable to express them in terms of the nodal DOFs (i.e., nodal displacements and nodal rotations ), which are gathered in the vector = [ 1 1 2 2 ⋯ N N ] T given as = . (7) is a square matrix, where each row with an odd number corresponds to the vector evaluated at node ( + 1)∕2, while each row with an even number contains the vector , evaluated at node ∕2. For the sake of clarity, the explicit expression for matrix is given here Solving the linear system of equations (7), we obtain the shape functions in terms of the nodal DOFs Finally, we substitute the solution for the coefficients -see Equation (9)-into Equation (5) and obtain the following relation which provides an explicit expression for the C 1 -continuous shape functions of an Euler-Bernoulli beam element The entity is commonly referred to as the matrix of shape functions. It contains 2 N shape functions of order = 2 N − 1. The numerical properties of these shape functions are highly dependent on the position of the nodes and therefore, we will discuss two different sets of nodes that are often used in FEA. In particular, equidistant (EQ) and GLL interpolation points are of interest. The latter choice leads to the C 1 -continuous SBEs as proposed in Ref. [74].

EQ points
In commercial finite element software, often only EQ (evenly-spaced) nodes are utilized and therefore, a brief description of this particular nodal distribution is called for. Note that in the remainder of this article, equidistant points are referred to as EQ points. Their coordinates can be easily computed by utilizing the expression below The shape functions for a 8-node element based on EQ points are depicted in Figure 2, which already demonstrate one important drawback of using this particular nodal distribution for interpolation purposes, that is, oscillations at the interval boundaries. These oscillations increase dramatically in magnitude with the number of interpolation points and lead to Runge's phenomenon as discussed in Section 3.1.

F I G U R E 2
Shape functions of a 8-node beam element based on EQ interpolation points.

GLL points
In the wide body of literature, it is generally acknowledged that the nodes should be clustered at the boundaries of the interpolation interval. This behavior is commonly ensured by selecting the roots of an appropriate family of orthogonal polynomials [56,88]. In the context of SEs, GLL and GLC 2 points are commonly employed. The coordinates of GLL points are defined as Here,̂denote the roots of the Lobatto polynomials or by exploiting Rodriguez' formula Based on the definition of the Legendre polynomials it is easy to obtain their first derivatives corresponding to the Lobatto polynomials −1 The sought after interpolation pointŝare exactly the roots of the polynomials defined by Equation (16). Hence, the complete set of interpolation points GLL is obtained by computing the zeros of the completed Lobatto polynomial of order + 1 defined aŝ+ 2 Gauß-Lobatto-Chebychev points: In numerical analysis, the properties of GLL and GLC points are virtually identical in most aspects. However, an accurate quadrature rule is only available for GLL points, while integration techniques relying on GLC points lack accuracy. Therefore, only GLL points are discussed in the article at hand.

F I G U R E 3
Shape functions of a 8-node spectral beam element based on GLL interpolation points.
The shape functions for an 8-node element based on GLL points are depicted in Figure 3. The behavior of the derived spectral C 1 -continuous shape functions, exploiting an Hermitian interpolation technique and GLL points, is distinctly different to that observed in Section 2.1.1. In this case, no oscillations are present and the maximum of the displacement shape functions is strictly 1 and occurs at its corresponding node. The ramifications of these observations are discussed in more detail in Section 3.

Semi-discrete equation of motion
As mentioned before, the derived SBE is based on the Euler-Bernoulli beam theory (considering only straight beam segments). Note that the beam-axis coincides with the global -axis. The material and geometric properties of the beam are given by Young's modulus , the mass density , the second moment of area of the cross-section , the cross-sectional area , and the length of the structure . In the Euler-Bernoulli beam theory, the axial displacement ( , ) approximated as where ( ) represents the (transverse) displacement of the beam. The axial (longitudinal) strain is related to the first derivative of the axial displacement ( , ) and consequently, given as The constitutive equations of linear elasticity for one-dimensional problems and a uni-axial normal stress state are exploited to compute the longitudinal stress as Equation (20) is derived by assuming a plane stress state and a negligible transverse normal stress . Using Equation (20), the bending moment bx can be related to the deflection ( ) as The variational formulation for transient beam problems is now obtained by utilizing Hamilton's principle. Using a standard Bubnov-Galerkin approach the semi-discrete equations of motion, corresponding to the weak form of the governing equations, read̈G where , , and denote the global (assembled) mass, damping, and stiffness matrices, G represents the vector of DOFs (deflections and rotations), and G is the load vector. The global quantities are assembled in a standard finite element procedure from the elemental ones, which are defined as wherêis the local strain-displacement matrix, defined according to Equation (19) aŝ and ( ( )) denotes the transverse normal force per unit length. Note that in all equations given above, the argument ( ) highlights the fact that the corresponding quantity depends on the global axial coordinate , which can be expressed in terms of the local coordinate for quadrature purposes. We can infer from the equations above that we do not only need the shape functions themselves, but also the second derivative with respect to the local coordinate . In analogy to Equation (11), the second derivative is defined aŝ( with ( ) , given by where each element of this vector can be compute by the following expression The definition of the damping matrix e is, at least in terms of a concise theoretical basis, not trivial at all and therefore, briefly discussed in Appendix A. Since there is no elegant mathematical derivation of the elemental damping matrix, a different (more phenomenological) approach has to be taken, which is based on Rayleigh's hypothesis. However, at this point we want to concentrate on mass lumping as an important means to realize efficient algorithms for explicit dynamics.

Mass lumping technique
Mass lumping is one of the most important topics in explicit dynamics, since it is a prerequisite for approaches based on the element-by-element (EBE) paradigm. Note that these types of approaches are essential for devising efficient (explicit) time stepping algorithms. In general, the semidiscrete equations of motion (22) can be re-written for many time stepping schemes in a generic, pseudo-static, way aŝ+ wherêis the effective or dynamic stiffness matrix and denotes the effective load vector. The superscript □ +Δ expresses that the corresponding quantity is evaluated at the discrete time + Δ , with Δ representing the selected time step size. Considering the central difference method (CDM), the effective quantities take the following form For more information on how to arrive at the above equations, please consult the monographs by Chopra [90] or Bathe [91]. Note that these two references are also an excellent source for a more detailed discussion on time stepping and the solution of transient problems in general. However, the important thing to realize from Equation (31) is that the dynamic stiffness matrix is only a function of and , but not (as would be the case for implicit methods). Consequently, the solution of the system of equations (30) is trivial, if both and are available in diagonal form. Therefore, mass lumping schemes play a critical role in explicit dynamics. Exploiting these advantages, a massively parallel implementation of the CDM is reported in Ref. [92], where examples with up to one billion DOFs are solved. Additionally, a thorough discussion of different mass lumping techniques as the enabling technology for explicit dynamics can be found in Refs. [93][94][95], where its application to Serendipity and spectral elements is discussed in detail.

2.3.1
Nodal quadrature The nodal quadrature technique is an important mass lumping scheme that is often also referred to as optimal lumping [96,97], since it is possible to obtain the theoretically optimal ROCs under certain conditions, as reported in Ref. [95]. Most discussions are, however, only related to C 0 -continuous continuum elements based on tensor product formulations of finite elements. Considering beam and shell elements the situation is quite different due to the existence of rotational DOFs. In the following, the basic idea of the nodal quadrature technique is recalled and its application to SBEs is discussed.
In cases, where the numerical integration of the mass matrix is accurate enough, the nodal quadrature technique is considered to be a variationally consistent approach to mass lumping, that is, convergence of the solution is still guaranteed, and therefore, this method is mathematically most appealing. In this method, the quadrature rule and the distribution of the nodes used to derive the shape functions are closely connected. To be specific, the basic idea consists in selecting an identical set of points for both the definition of the quadrature rule and the construction of the polynomial interpolant, serving as the basis/shape function. Generally speaking, mass lumping by nodal quadrature fulfills the conservation of mass requirement, but cannot ensure the positive-definiteness of the mass matrix [98][99][100]. This is no problem for continuum SEs [95], but constitutes a severe problem for its application to C 1 -continuous beam and shell element formulations. Moreover, it is known that the optimal order of accuracy is only achieved in the numerical analysis, if the applied ( + 1)point nodal quadrature rule is exact for polynomials of order 2 − 1 or higher [101]. That is to say, a slight under-integration of the highest order terms in the mass matrix integrand, which are of order 2 , is allowed. This condition is fulfilled by the GLL quadrature rule, if + 1 integration points are used per direction 3 , which is generally the case for continuum SEs. Despite the fact that the maximal polynomial order of the integrands to compute the elemental mass and stiffness matrices for SBEs-Equations (23) and (24)-is still 2 , the number of integration points is reduced to only ( + 1)∕2 (recall that we have to fulfill two conditions at each node to ensure C 1 -continuity of the shape functions), which leads to a severe under-integration, which is in stark contrast to C 0 -continuous continuum SEs (see Table 1). This error might be nonetheless acceptable for the mass matrix of SBEs, since we gain an increased efficiency by having a diagonal formulation readily available. Considering the computation of the stiffness matrix e and the load vector e , there are, however, no additional advantages to be gained and therefore, a conventional Gaussian quadrature scheme-( + 1)-point rule-is employed 4 .
By definition, the spectral beam shape functions hold the Kronecker delta and partition of unity properties. These features are indispensable for the nodal quadrature method to work. It means that all off-diagonal components of the mass matrix are zero and only the diagonal components exhibit values that are different from zero. However, not only the off-diagonal components are zero but also those diagonal components associated with the rotational DOFs (see Figures 2b 3 The order of accuracy of a GLL quadrature rule is 2 N − 3 and for C 0 -continuous shape functions the polynomial degree is given by = N − 1. Consequently, the order of accuracy of a GLL quadrature rule in terms of is 2 − 1. 4 Considering the load vector, a (( + 1)∕2 + 1)-point rule would be theoretically sufficient, assuming that the distributed load ( ) is uniform along the length of the SBE. TA B L E 1 Relation between the number of nodes N , polynomial degrees of the shape functions SF and the integrand M to compute the mass matrix for spectral continuum and beam elements, and the order of the polynomial int that can be exactly integrated by a N -GLL rule a .

Continuum
Beam Continuum Beam   2  1  3  2  6  1   3  2  5  4  10  3   4  3  7  6 1 4 5 5 4 9 8 19 7 a Remark: For continuum elements (C 0 -continuous shape functions), the relation between the number of nodes N and the polynomial order of the shape functions is: SF = N − 1. Considering beam elements (C 1 -continuous shape functions), this relation is different and can be expressed as: SF = 2 N − 1. This is due to the fact that we have to account for both displacement and slope continuity in the derivation of C 1 -continuous shape functions. However, we only have N integration points at our disposal independent of the selected element type, since the nodal quadrature technique is implemented to achieve a diagonal mass matrix. Using this approach, a polynomial of order int = 2 N − 3 is integrated exactly. In terms of the polynomial degree of the shape functions, we notice that int = 2 SF − 1 for continuum elements. In contrast, for beam elements int = SF − 2 holds. This is one critical aspect in the application of the nodal quadrature technique for SBEs. and 3b). This renders the lumped mass matrix positive semi-definite, which is unfortunately not compatible with explicit time integration schemes. Before we discuss possible remedies in Section 2.3.2, let us recall the expression for the lumped mass matrix. In terms of numerical integration, the elemental mass matrix, given in Equation (23), can be re-written as

Treatment of zero masses
As mentioned before, it has been found that the lumped mass matrix exhibits zero components associated with the rotational DOFs. In Ref. [74], two methods to circumvent this drawback have been proposed. For the sake of completeness, these approaches are recalled at this point and will be thoroughly assessed in Section 3 regarding their influence on the quality of the solution compared to simulations with a consistent mass matrix.

Discrete mass method
The discrete mass method (DMM), in Ref. [74] referred to as direct explicit method (DEM), is a heuristic approach, where zero diagonal entries in the mass matrix are replaced by a small value 0 , defined as The scaling parameter is a large constant and a value of 10 7 has been suggested in previous studies. Using this approach the diagonal mass matrix becomes positive-definite and can be directly employed in conjunction with Equation (30) to advance in time.

Static condensation
The static condensation (SC), which is a special case of the Guyan reduction [102], is a dimensionality reduction method, which reduces the overall number of DOFs. This is achieved by separating the DOFs into master (m) and slave (s) sets 5 .
The slave DOFs are condensed in the reduction process, which leaves us with a system of equations only containing the master DOFs. Thus, Equation (22) In the next step, we can enforce the equation in the second row of Equation (35) for s G as Equation (36), as given above, is the exact expression for arbitrary numerical methods. However, by exploiting the particular structure of the mass matrix of the proposed SBEs, Equation (36) can be further simplified. For any diagonal mass matrix formulation, we know that ms = sm = holds. Moreover, all slave DOFs are also associated with zero mass components for the developed SBEs, meaning ss = , and therefore, we obtain Note that so far no approximation has been introduced with respect to the proposed SBEs and therefore, no additional errors are caused by this procedure. This result is substituted into the first row of Equation (35) mm̈m and we obtain mm̈m As we can clearly observe, Equation (35) has been successfully re-written only in terms of the master DOFs only. In a more compact notation, Equation (39) is given as mm̈m G +̂m m m G =̂m G , witĥm m = mm − ms ( ss ) −1 sm and̂m G = m G − ms ( ss ) −1 s G .
Equation (40) can be directly used in conjunction with the CDM to compute the time history response of the master DOFs, while the values for the slave DOFs are determined in a post-processing step utilizing Equation (37). Note that if no loads are acting on the slave DOFs,̂m G can be simply replaced by m G . At this point, we want to highlight again that the introduced static condensation is not approximate in nature (when suing the proposed SBEs), since a lumped mass matrix is readily available and additionally, the slave DOFs do not contribute any mass to the system. Consequently, all errors that are observed when applying this approach are directly related to the fact how well the lumped mass matrix is capable of representing the consistent one. In this context, the static condensation can be seen as an exact method, whereas the DMM is only approximate in nature and the best results it can deliver are those obtained by the static condensation.
For the sake of completeness, we want to provide the complete expression for reducing the number of DOFs in transient analysis for general numerical methods, where only a consistent mass matrix is available. Additionally, we will include physical damping in the equations of motion 6 . One of the main assumptions in the static condensation approach with respect to general numerical methods, is that the inertia terms are negligible and that no loads are applied at the slave DOFs. In this case, the transformation matrix , reducing the full system of equations can be given as which is equivalent to Equation (37). Thus, the reduced system of equations can also be written in an alternative form as 6 Remark: Strictly speaking, none of the two approaches for treating zero masses must be utilized for the proposed SBE, since the damping matrix effectively regularizes the dynamic stiffness matrix̂-see Equation (31). Thus, the selected time stepping method can be directly applied to the derived system matrices.
Finally, we remind the reader of the fact that the introduced reduction technique does not introduce errors in static analyses. The same comment applies for numerical methods featuring a diagonal (lumped) mass matrix and zero masses at the slave DOFs, which is the case for the proposed SBEs, when no physical damping is present. In other words, the static condensation is exact for the element type investigated in this contribution both in static and (undamped) dynamic analyses. However, for other numerical methods the underlying assumptions generally cause deteriorate results in the high-frequency range as only frequencies close to the lowest eigenfrequencies of the system can be expected to be accurately captured.

NUMERICAL PROPERTIES
In this section, important numerical properties of the high-order SBEs, such as interpolation, conditioning, and numerical ROCs are investigated in detail. The accumulated findings will guide the application of the novel element type and provide us with valuable insights regarding advantages and possible limitations of this formulation.

Interpolation properties and Runge's phenomenon
Ideally, High-order polynomial shape functions should be capable of interpolating smooth solution fields to an arbitrary degree of accuracy, which means the error should decay until machine precision is reached. Thus, the interpolation properties of shape functions are an important aspect that needs to be considered in order to assess their performance. To illustrate these properties, Runge's function which is a smooth function within the interval [+1, −1], is often employed (see Figure 4). Using any available interpolation scheme, it is usually expected that the quality of the approximation of the original function is improving with an increased number of interpolation points. However, the distribution of these interpolation points cannot be chosen arbitrarily, as will be shown in the following for both Hermitian and Lagrangian interpolations.
In the classical FEM, it is customary to use EQ points to develop high-order elements, while GLL-or GLC-points are favored in the SEM. Generally, it is expected that an increasing number of interpolation points invariably leads to an improved quality of the approximation of the original function. The opposite is, however, the case for EQ points, where the magnitude of the oscillations near the boundaries of the interpolation domain is increasing with the number of interpolation points (see Figure 5), which will eventually cause the interpolation to fail. This effect clearly illustrates  that choosing an evenly-spaced set of interpolation points is not only not optimal to achieve high-quality results, but it is a fundamentally wrong approach. The described behavior can be easily circumvented by employing non-equidistant points clustered at the interval boundaries [56], such as GLL-and GLC-points. Since the overall properties of both sets of interpolation points are nearly identical, we restrict our discussion in the remainder of this article to GLL-points. Note that such distributions decrease the Lebesgue constant, which is a measure of how suitable the interpolant of a function (at the chosen nodal distribution) is in comparison with the best polynomial approximation of the function. In this sense, the Lebesgue constant also bounds the interpolation error.
In Figure 5, the interpolation results for Runge's function are depicted. As mentioned before, it is observed that EQ points inevitably lead to oscillations, while this problem is suppressed when utilizing GLL points. As the number of interpolation points is successively increased, a monotonous convergence of the error is seen in the case of GLL points. In contrast, interpolations based on EQ points do not show a monotonic convergence at all. For our assessment, the error is calculated in terms of the area below the interpolated and original curves as When employing GLL-points spectral convergence is guaranteed, which means that the interpolation error decreases with the number of interpolation points as In comparison to the Lagrangian interpolation, we note that the Hermitian interpolation based on GLL-points is better behaved, meaning that it results in a reduction of the error for the same number of interpolation points, while also minimizing still present oscillations. Overall, it can be stated that the proposed Hermitian interpolation technique achieves similar results to the Lagrangian interpolation scheme. The latter approach has been thoroughly studied in the context of SEM and at this point, we conjecture that the positive properties noted for Lagrangian interpolation schemes carry over to the new method.

Condition number of the stiffness and mass matrices
The conditioning of the system matrices is of utmost importance for the solution of linear systems of equations or eigenvalue problems that arise in static, transient, or modal analyses. It is well-known that the condition number is one of the main factors influencing the number of significant digits for direct solvers or the required iteration count for iterative ones. Thus, we need to investigate how the condition number evolves with an increasing order of the utilized shape functions. To this end, condition number estimates known from C 0 -continuous shape functions are taken and adapted to the Hermitian interpolation employed in this research. A comprehensive discussion on the conditioning of the stiffness matrix obtained from different numerical methods, such as SEM, p-FEM, and IGA, can be found in Ref. [103] and the references cited therein. The numerical results, which are discussed in the remainder of this section, confirm again that GLL points show distinct advantages over EQ points, and are thus generally preferable also in Hermitian interpolation techniques. In the context of this paper, we are only interested in the evolution of with respect to a p-extension, that is, the polynomial order of the shape functions is elevated, while the element diameter (size) ℎ is fixed. On the other hand, for h-extensions, where the polynomial order of the shape functions is fixed, while the element diameter is decreased, it is well-known that the condition number scales with ℎ −2 for a wide variety of element types [104,105]. Note that this estimate holds independent of the dimensionality of the problem. Unfortunately, the condition number is much more sensitive to the polynomial degree of the shape functions, as will be discussed in the remainder of this section. Considering C 0 -continuous SEs based on Lagrangian shape functions defined on GLL points, Maitre and Pourquier, Melenk, and Gervasio et al. derive the following estimate [106][107][108] This estimate is shown to be sharp not only for single undistorted elements but also for general discretizations [103], which contradicts the conjecture provided by Olsen and Douglas Jr. in Ref. [104], stating that p-version matrices are bounded from below by 4 for any dimensionality. However, their estimate for the spectral condition numbers of high-order shape functions of Lagrange-type defined on equidistantly spaced points being is very accurate [104]. The expressions provided in Equations (46) and (47) serve as a basis for assessing the conditioning of the proposed Hermitian SBEs. These formulas are taken as starting point and it is conjectured that the condition numbers for the proposed element behave in a similar way. Consequently, a generic formula describing the relation between condition number and polynomial degree takes the following form where, the introduced constant includes the dependency on the element size ℎ, while might be a function of the dimensionality of the problem , if an extension to plate and shell elements is contemplated. Note that the formulation of two-and three-dimensional SBEs requires only the multiplication of the element matrices with a suitable transformation matrix, which does not have an influence on the condition number. In this sense, refers to the dimensionality of the element in its local coordinate system. Consequently, the analysis of one-dimensional examples is sufficient to study the numerical behavior.
To determine whether Equations (46) and (47) are suitable for describing the condition number of elements based on Hermitian shape functions defined on GLL or EQ points, respectively, we can compute a best fit of our numerically obtained data. To this end, the data based on GLL points is fitted to while a different expression is used for the data based on EQ points log 10 ( EQ ) − log 10 (4) = (ℎ) + ( ) log 10 ( ).
Comparing Equations (49) and (50) to Equations (46) and (47), we can identify the values for and when continuum SEs are employed. When GLL points are selected, (ℎ) = ℎ −2 and = 3. On the other hand, choosing EQ points leads to the same (ℎ), but a different exponent = 2 − 4 . Note that the term 4 is dominating the value of the condition number for EQ points. In our fitting procedure, this term has been moved to the left-hand side of Equation (47). So far, only condition number estimates for the stiffness matrix have been reported. In Ref. [108], also an estimate for the condition number of the mass matrix based on GLL points is provided as Due to the similarity of the estimates for GLL points with regard to the condition number of the stiffness and mass matrices, and the fact that there is no estimate (at least to the authors' knowledge) for the conditioning of the mass matrix based on EQ points, the mass matrices will also be fitted to the expression provided in Equations (49) and (50).
In the remainder of this section, we will only report the values for as the dominating factor for the numerical behavior of the condition number. The constant is of less interest as it only causes a parallel shift of the curve. In order to assess whether the estimates derived for C 0 -continuous shape functions also hold for the proposed C 1 -continuous Hermitian shape functions, a simply-supported beam (i.e., the transverse displacement is set to zero at both ends of the beam) is discretized by a single SBE and the number of nodes N is successively increased from 2 to 10, with = 2 N − 1. For all discretizations, the condition numbers with respect to the stiffness matrix ( ), the mass matrix ( ), and the product −1 ( −1 ) are computed. The conditioning of the last term is especially important for the stability limit in explicit dynamic analyses.
The numerical results with respect to the conditioning of the elemental matrices are depicted in Figure 6, where both the actual recorded data (red diamond-shaped markers: ◊) and the best fit curves (solid black line: ) are included. The values reported for the mass matrix are related to the consistent formulation, that is, a full numerical integration, where a standard Gaussian quadrature rule is employed. This is done for a more meaningful comparison since no accurate nodal quadrature rule is available for EQ points. We observe a good agreement between the condition numbers of the stiffness F I G U R E 6 Condition numbers of the system matrices ( , , and − ) versus polynomial order of the C 1 -continuous Hermitian shape functions based on EQ or GLL points for a simply-supported spectral beam element.

F I G U R E 7
Condition number of the matrix − versus polynomial order of the C 1 -continuous Hermitian shape functions based on EQ points for a simply-supported spectral beam element-the fitting is based on Equation (49) instead of Equation (50). matrices for both nodal distributions (see Figures 6a and 6b) and the estimates proposed in Equations (49) and (50). The obtained -values are 6.449 and −1.012 for the GLL and EQ points, respectively. For one-dimensional continuum SEs, the theoretical -values are 3 and −2 (see to the results published in Ref. [103]). These results confirm that the conditioning of SBEs is worse compared to that of spectral bar elements.
The fitting results for the condition number of the mass matrix are depicted in Figures 6c and 6d. The agreement is still quite good for the mass matrix based on GLL points, but for EQ points a different estimate needs to be derived, which is however, out of the scope of this contribution. For the GLL case, the -value is determined as 5.902, which is again significantly higher compared to the theoretical value of 4 given in estimate Equation (51). To summarize, the results obtained so far, highlight again that using EQ points leads to a significant increase (by several orders of magnitude) in the condition number of the system matrices.
In Figures 6e and 6f, the fitted curves for the matrix product − are shown. The numerical behavior of the condition number is very similar to what has been noted for the mass matrix. Considering GLL points, an -value of 6.560 is obtained. By closer inspection, we note that the curves for GLL and EQ points are very similar and therefore, the conditioning data for the matrix product − based on EQ points is fitted again using Equation (49) (see Figure 7, -value: 6.556). This results in a much better agreement compared to the curve shown in Figure 6e, where Equation (50) has been employed. Surprisingly, it must be stated that the conditioning for − seems to be independent of the choice of interpolation points. At this point, the reason for this behavior is not yet understood and requires additional investigations.

Numerical rates of convergence
The numerically attainable ROCs are a good indicator of the overall performance of a novel element technology. Therefore, we will study the convergence of the proposed SBEs considering static, modal, and transient analyses. In the remainder of this section, only academic benchmark examples are investigated, while more complex problems are tackled in Section 4.

Static analysis
In this section, we evaluate the performance of the proposed SBE for static applications to demonstrate that the theoretically optimal ROCs are achieved. To this end, we will investigate a simple model problem, which consists of a cantilever beam that is subjected to a distributed transverse load as depicted in Figure 8. Although it is easily possible to derive an analytical solution for this particular example, it is nonetheless worth studying, as we will observe many of the features which also appear in more complex models. Especially, the concept of p-refinement and its consequences for the quality of the solution in comparison to an h-refinement are readily illustrated in a one-dimensional setting. First, let us recall the fourth-order ordinary differential equation governing beam bending problems This formulation is particularly helpful to analytically determine the (transverse) displacement ( ) of a beam structure subjected to an arbitrarily distributed normal load ( ) and any number of point forces . Note that Equation (52) is straightforwardly derived from Equation (21) by double differentiation.
As mentioned before, we employ a cantilever beam, which is subjected to a distributed load without applying any point forces (see Figure 8 for a sketch of the model). The load is given as a trigonometric function ( ) = sin( 1 ) to ensure that the analytical solution is not part of the ansatz space of the SBEs. Hence, the closed form solution to Equation (52) is given as where the values of the integration constants are related to the specific BCs. The derived analytical solution is based on the assumptions that the material is linear elastic and homogeneous, and that the beam features a constant cross-section, that is, and are no functions of the coordinate . For a cantilever beam 7 , the following BCs are applicable with 1 being the coordinate of the left end of the beam and 2 denoting the coordinate of the right end. Hence, the length of the beam is given as Consequently, the exact solution of our simple model problem reads In order to assess the quality of the numerical solution, we compute the 2 -error in the displacement field as 7 Remark: The exact dimensions of the beam and its specific cross-section do not influence the convergence behavior demonstrated in this section and therefore, we refrain from providing exact values for the beam parameters, unless explicitly stated otherwise.

F I G U R E 9
Convergence curves for the spectral beam element in static analysis-p-and h-extensions.
where the subscripts □ ex and □ num denote the exact (reference) and numerically obtained solutions, respectively. Note that the distributed load ( ) results in a smooth displacement field-Equation (57)-and therefore, optimal convergence of the solution is expected. The convergence plots for the static example are depicted in Figure 9, for both EQ and GLL points. We observe that the expected exponential convergence is achieved when a p-extension is applied, that is, the number of nodes N is successively elevated from 2 to 9. In this case, the structure is only discretized by a single SBE. On the other hand, an algebraic convergence with optimal slope is obtained for h-refinement approaches, where the number of elements is doubled for each simulation starting with 1 element. For example, a 2-node SBE exhibits a decay of the error with (ℎ −3 ), where ℎ denotes the element size. Optimal ROCs are also attained for a 5-node SBE with (ℎ −9 ). Note that the theoretically predicted ROCs are (ℎ − ). Apart from worsening results for the 9-node SBE when EQ points are employed, we see an excellent agreement in the results independent of the selected nodal distribution. Therefore, we can conclude that in static applications, when using a moderate number of nodes, the actual nodal distribution does not exert a significant influence on the quality of the results. It must be stressed that a different behavior will be observed for elements featuring a larger number of nodes, where Runge oscillations deteriorate the approximation/interpolation properties of the shape functions based on EQ points. In numerical examples, this behavior is hard to show since already 5-node SBEs yield very accurate results for all problems exhibiting smooth solution fields.

Modal analysis
In this section, the eigenvalues and mode shapes of selected structures are computed, which is commonly referred to as modal analysis. The modal analysis problems, solved in the following, will give us a first insight into the influence of the mass lumping procedure on the accuracy of the solution. Here, it is hoped that the lumped mass matrix formulation does not result in (significantly) deteriorated results compared to the consistent (fully-populated) mass matrix. To this end, the eigenvalues of simply-supported and cantilever beams (see Figure 10) are computed numerically and compared with readily-available analytical solutions. As is well-known, the eigenvalues are related to the eigenfrequencies of the system F I G U R E 1 0 Sketch of the beam structures subjected to the modal analysis. by = 2 . In this section, the mass matrix is lumped by nodal quadrature, as described in Section 2.3.1, but the zero masses are not explicitly treated (see Section 2.3.2), which results in infinite eigenvalues. These solutions are consequently discarded in the subsequent analysis of the element's performance. The error in the eigenvalues can either be determined as the relative error in the EV -th eigenvalue = | EV ,ex − EV ,num | EV ,ex (59) or by an average error over the first EV eigenvalues For the current assessment, the first 20 eigenvalues ( EV = 20) are taken into consideration. Since the overall error-value is usually dominated by the higher eigenvalues, as they are less accurately resolved by the spatial discretization, we rely on Equation (59) in the assessment of the numerical accuracy. It has been noted that the obtained convergence is very similar when employing Equation (60). In the next two paragraphs, the analytical solutions for the two suggested models are recalled before the numerical results are presented.

Simply-supported beam
For a simply-supported beam with constant cross-section and homogeneous material properties, the circular eigenfrequencies are given by while the mode shapes follow a simple trigonometric function where is an arbitrary constant which can be used for normalization purposes. As an example, the first 20 mode shapes have been depicted in Figure 11.

Cantilever beam
For a cantilever beam with constant cross-section and homogeneous material properties, the circular eigenfrequencies are principally also given by Equation (61) but the parameters are defined differently and must be determined by solving the following transcendental equation For the first 20 eigenfrequencies, the values for are listed in Table 2. Compared to the simply-supported case, also the mode shapes are more complex for the cantilever beam ] .
Due to the fact that the solutions for the respective mode shapes-the first 20 are visualized in Figure 12 are again smooth functions, an optimal convergence of the results is expected, if a consistent mass matrix (fully-integrated and -populated) formulation is employed. Therefore, we can easily study and observe the effect of mass lumping by nodal quadrature (under-integration) on the eigenvalues, since there are no other effects present that could deteriorate the numerical results.

Numerical results
The results of the numerical analysis of the attainable ROCs are depicted in Figures 13 and 14 for the simply-supported and cantilever beams, respectively. The results are related to the consistent mass matrix implementation of the beam element.
For both example problems, we observe the predicted exponential convergence when conducting a p-refinement, where the number of elements is selected as E = 10, while the number of nodes is elevated from 2 to 9. An additional h-extension with 10, 20, 40, and 80 elements reveals that the theoretically optimal ROCs of order (ℎ −2 ) are recovered initially for both example problems and nodal distributions. Overall, the agreement is better for the cantilever beam example since the solution is more complex. In the case of a consistent mass matrix, we have to state that despite a difference in the minimum error that is reached, we cannot report a significant difference in the performance of beam elements based on EQ or GLL points. The reader has to keep in mind, however, that this statement is related to a moderate number of nodes (two to five nodes) only, where Runge's phenomenon, discussed in Section 3.1, is not very pronounced. On the other hand, a 5-node SBE already features shape functions of order 9, which is sufficiently high for most applications of practical interest.
Considering an implementation of the SBEs in conjunction with a lumped mass matrix formulation, where the zero masses are treated by the DMM (see Section 2.3.2), the picture is quite different. Note that although it is theoretically also possible to diagonalize the mass matrix of SBEs based on EQ points by means of Newton-Cotes quadrature rules, we observe zero and negative masses and therefore, refrain from applying this technique. This approximation of the mass matrix invariably yields diverging results in transient analyses and consequently, we only discuss the results obtained by the SBEs based on GLL points and GLL quadrature rules. The results for both example problems are reported in Figure 15 F I G U R E 1 2 Mode shapes of the clamped-free beam.

F I G U R E 1 3 Convergence curves for the spectral beam element in conjunction with a consistent mass matrix formulation in modal analysis (simply-supported beam)-p-and h-extensions.
for a -value of 0, that is, the rotational masses are not treated and are kept as zero-see Equation (34). Although, we observe a convergence of the results, we also have to admit that the results are worse than initially expected. The attainable error is sufficient from an engineering point of view, but orders of magnitude higher compared to the consistent mass matrix case. The deteriorate results are solely attributed to the severe under-integration of the mass matrix, due to the fact that no additional approximations have been included in the formulation. That is to say, the lumped mass matrix (obtained by the nodal quadrature technique) is not capable of approximating the consistent mass matrix to a high degree of accuracy. However, on the plus side, we at least observe a convergence of the numerical eigenfrequencies to the analytical values, which is not possible to obtain when using EQ points.

F I G U R E 1 4 Convergence curves for the spectral beam element in conjunction with a consistent mass matrix formulation in modal analysis (cantilever beam)-p-and h-extensions.
F I G U R E 1 5 Convergence curves for the spectral beam element (GLL nodes) in conjunction with a lumped mass matrix formulation and the use of the DMM for treating zero masses in modal analysis-10 spectral beam elements, p-extension.
Considering transient problems, it is not possible to use a value of 0 for in combination with explicit time integration schemes. Therefore, Figure 16 depicts the convergence behavior for a selection of different -values in the modal analysis of the cantilever beam. It is observed that for > 1 × 10 5 a good agreement of the results is achieved. Only for the case of 2-node beam elements, we notice an increased error that is directly related to the value of the added discrete mass. This effect is associated with the increasing difference in magnitude of the masses corresponding to the displacement and rotation degrees of freedom, respectively. The higher the -value, the smaller the added rotational mass. It is well-known that such a difference in matrix components will ultimately lead to a conditioning problem. Increasing the value of beyond 1 × 10 7 , does not seem to have any further positive effects on the attainable accuracy and therefore, = 1 × 10 −7 is chosen for all further examples.
Instead of employing the DMM to treat the zero masses in the lumped mass matrix, we can also employ the Guyan reduction technique (see Section 2.3.2) to condensate all rotational DOFs from the system. The advantage of this approach is that no additional parameter has to be chosen. Note that for the proposed SBEs, the Guyan reduction approach is exact for both static and dynamic problems, due to the properties of the lumped mass matrix. This is in stark contrast to other F I G U R E 1 6 Influence of the -value (DMM) on the numerical convergence for the spectral beam element (GLL nodes) in conjunction with a lumped mass matrix formulation and the use of the DMM for treating zero masses in modal analysis (cantilever beam example)-p-extension. numerical methods, where the condensation is only exact for static analyses, but introduces additional errors in dynamic analyses. Therefore, we must stress again that all errors are related to the lumping of the mass matrix by under-integration. In Figure 17, the convergence curves for the two examples are shown. We observe that the convergence curves are in virtually coincident with those obtained for the DMM case. That is to say, independent of the method for treating the zero masses, a convergence of the numerical results to the analytical solution is observed. However, again we must stress that only sub-optimal ROCs are obtained, which are far from the exponential convergence noted for the consistent mass matrix formulation of the SBEs. However, it is a rather interesting finding and confirms that the numerical results for both techniques to treat the zero rotational masses are essentially identical. At this point, we can conjecture that the use of either the DMM or Guyan reduction is equivalent when conducting modal analyses.
To get a better picture of the convergence behavior for the lumped mass matrix case, the number of elements is increased to E = 50. This spatial discretization is fine enough to reveal the asymptotic ROC. In Figure 18, the convergence curves when using the DMM are shown, while Figure 19 contains the plots for the Guyan reduction case. The convergence F I G U R E 1 7 Convergence curves for the spectral beam element (GLL nodes) in conjunction with a lumped mass matrix formulation and the use of the Guyan reduction for treating zero masses in modal analysis-10 spectral beam elements, p-extension.

F I G U R E 1 8
Convergence curves for the spectral beam element (GLL nodes) in conjunction with a lumped mass matrix formulation and the use of the DMM for treating zero masses in modal analysis-50 spectral beam elements, p-extension.

F I G U R E 1 9
Convergence curves for the spectral beam element (GLL nodes) in conjunction with a lumped mass matrix formulation and the use of the Guyan reduction for treating zero masses in modal analysis-50 spectral beam elements, p-extension.
properties are again essentially identical and therefore, our previous statement that in a modal analysis it does not matter which approach is employed is confirmed. In both cases, the error decays at an order of around (ℎ −5 ). Hence, we clearly do not obtain an exponential convergence, which is the case for the consistent mass matrix formulation. The reduction in the attainable ROCs is the price we have to pay for diagonalizing the mass matrix. So far, the use of the Guyan reduction (static condensation) technique is recommended due to the fact that no additional parameter has to be carefully selected.

Frequency spectrum of a simply-supported beam
In the following, the discrete spectrum is computed for the simply-supported beam problem (described above) with different spatial discretizations. This approach is analogous to the one proposed in Ref. [109] in the context of IGA. To this end, the number of nodes per element is increased from 2 to 5, while the number of elements is adjusted such that models with 1,000 DOFs are created. Note that for the case N = 4, a numerical model featuring 1,002 DOFs and 167 spectral beam elements is generated. For all models, the entire spectrum, that is, 1000 modes are computed and compared to the analytical solution given by Equation (61). The obtained results are presented in Figure 20, where the numerical eigenfrequencies which is plotted over the normalized mode number rel , which is obtained by dividing the by the number of DOFs of the structure under consideration. We observe that also for C 1 -continuous elements based on Hermitian interpolation schemes both acoustical and optical branches exist. Depending on the number of nodes of the SBE, N − 1 optical branches exist in the discrete spectrum if a consistent mass matrix formulation is employed. This behavior is very similar to standard C 0 -continuous displacement-based finite elements, where the number of optical branches is also related to the number of nodes in an element, that is, N optical branches exist in the discrete spectrum. Consequently, less optical branches are observed when employing the proposed SBEs. At this point, we have to emphasize that optical branches in the discrete frequency spectrum of the structure are directly related to spurious waves and therefore, need to be either physically or numerically damped out of the response of the structure if excited. On the other hand, when using a lumped mass matrix formulation based on the nodal quadrature technique the results are much worse, especially in the high-frequency regime. Due to the under-integration, which is introduced to achieve a diagonal mass matrix, a higher frequency error has to be accepted and also the behavior is more erratic and less predictable. While an error of well below 5% can be easily achieved by the consistent mass matrix formulation for half the spectrum, an error of above 5% is already reached at a quarter of the spectrum for the lumped mass matrix case with elements having more than two nodes. To be more precise, the error threshold is first exceeded at the following values for rel : Consequently, one has to be very careful that the incorrectly approximated frequencies are not part of the solution spectrum when lumping is required for the sake of efficiency of the explicit solution procedure. Due to the fact that isogeometric elements are capable of providing a higher inter-element continuity (up to C −1 ) depending on the polynomial degree of the shape functions [109], no optical branches are encountered 8 , which is in contrast to the proposed SBEs. Nevertheless, the observed behavior is still an improvement compared to classical C 0 -continuous displacement-based finite elements (see also discussions in Ref. [110]).

Transient analysis
In this section, a transient example is investigated to further assessment the numerical performance of the proposed SBEs. In particular, a wave propagation example with an analytical solution has been selected. Here, a simply-supported beam with a prescribed initial displacement field is considered [74]. The spatial variation of the initial displacement field is given as while its first spatial derivative (slope), which is needed to prescribe the rotational DOFs of the SBEs, reads The amplitude of the prescribed displacement field 0 and the parameter , controlling the width of the spatial pulse, can be chosen arbitrarily. As mentioned before, this problem can be solved in closed-form and the analytical solution is given in terms of a Fourier series as In Figure 21, the first four integrand functions (normalized to a maximum value of 1) to compute the analytical-see Equation (69)-are depicted. To obtain a reference solution for our numerical simulations, the Fourier series has to be truncated after a finite number of terms max . The number of terms has been chosen according to the value of̂, which is essentially the amplitude of the -th Fourier term. It is well-known, that the magnitude of̂is getting smaller with 8 Conjecture: The spectrum of C 1 -continuous isogeometric beam elements is identical to that of the proposed SBEs based on an Hermitian interpolation scheme, since an identical ansatz space is spanned.
increasing and therefore, the contribution of higher Fourier terms becomes negligible at some point. In this particular example, we decided to truncate the Fourier series, if the ratio of the largest to the smallest Fourier coefficient is below machine precision. Thus, the number of Fourier terms is set to max = 600. We also want to point out that̂= 0, when assumes an even number. Therefore, the values for ∈ {2, 4, 6, …} do not need to be computed explicitly, but they are directly set to 0. As can be inferred from Figure 21, the integrals for = 2 (red curve) and = 4 (curve) must be zero due to the point symmetry with respect to = 1 m.
To obtain the numerical solutions, presented in the following paragraphs of this section, the beam parameters (rectangular cross-section) are chosen as: length = 2 m, width = 100 mm, height ℎ = 2 mm, Young's modulus = 70 GPa, Poisson's ration = 0.3, and mass density = 2700 kg∕m 3 . In Figure 22, the beam model is sketched in its initial configuration, where the displacement field is magnified by a factor of 100. The initial spatial discretization consists of E = 40 SBEs, while the number of nodes is again successively elevated from 2 to 9. The results of this p-refinement process are depicted in Figure 23 for EQ and GLL points for a simulation time of sim = 500 s. Note that in the case of EQ points, results have been only obtained for the consistent mass matrix formulation, while also lumped mass matrix results are available for GLL points. The error has been measured in the 2 -norm for the displacement at the observation point , F I G U R E 2 2 Sketch of the beam structure subjected to the transient analysis: Simply-supported beam with initial displacement field (shown for = 2 m, 0 = 2 mm, and = 7.5 × 10 −3 1∕m). For the purpose of visualization, the displacements are scaled by a factor of 100.

F I G U R E 2 3
Convergence curves for the spectral beam element in transient analysis-p-extension.
First, we consider the results for both nodal distributions in conjunction with a consistent mass matrix formulationsee Figure 23. Again, we observe the expected exponential convergence with a slight advantage in terms of the attainable accuracy for the implementation based on GLL points. However, when it comes to the use of diagonal mass matrices, convergence is only achieved for GLL points. Overall, the numerical results show a similar trend as discussed in Section 3.3.2. As soon as a lumped mass matrix is employed, a sub-optimal convergence behavior has to be accepted, whereas it does not make any difference whether the DMM or Guyan reduction approach is taken. In Figure 24, the results of an additional h-refinement procedure for the models based on GLL points are depicted. For the consistent mass matrix case, the error decays with an order of (ℎ −( +2) ), which is actually higher than expected. On the other hand, the attainable ROC of the lumped mass matrix approach is always −5, independent of the number of nodes and thus, the polynomial order of the used shape functions. These results are only given for the Guyan reduction, since it has been shown so far that almost identical values are achieved by the DMM approach. Furthermore, we notice that an error plateau is reached at an error level of around 10 −5 %. Note that if sufficiently fine spatial discretizations are investigated, all convergence curves would level off at this error value. The reason is not related to the limited accuracy of the proposed SBEs, but is seen in the limited accuracy of the chosen time-stepping method. In all transient analyses, if not stated otherwise, we employ the trapezoidal rule form the Newmark family of time integration schemes, which is a second-order accurate scheme featuring no amplitude decay. To reach lower error levels, one would have to decrease the time increment from Δ = 1 × 10 −8 s to Δ = 1 × 10 −9 s or use high-order time stepping schemes, such as the one proposed by Song et al. [111]. This is, however, out of the scope of the present contribution and therefore, we only mention these possibilities without further pursuing them. Note that an implicit time integration scheme has been selected to accommodate the use of both consistent and lumped mass matrices. Moreover, it is not conducive in term of comparability of the obtained numerical results to employ different time stepping methods. This is by no means a limitation of our approach, since we are not interested in studying the accuracy of the temporal discretization. Readers interested in results obtained by a combination of the prosed SBEs and the central difference method are referred to the articles by Jain and Kapuria [74][75][76].
In Figure 25, the displacement field over the entire length of the beam is depicted for different time instances. Here, a relatively coarse spatial discretization has been chosen with 40 beam elements having 5 nodes. While the analytical solution (black curve) and the numerical results obtained with a consistent mass matrix formulation (red curve) are in excellent agreement (i.e., the error is well below 1%), we observe distinct differences when using a lumped mass matrix (the error is still above 20%). As discussed before, the DMM and Guyan reduction techniques provide similar results in terms of accuracy, but we note an increased occurrence of spurious oscillations when the DMM approach is utilized. This is in agreement with the fact that the Guyan reduction is an exact approach for the proposed SBEs, while the DMM technique is only approximate in nature. Therefore, differences are expected for low-resolution meshes. In Figure 26, the results for a finer spatial discretization featuring 80 beam elements having 8 nodes are depicted. In this case, the error is well below 1% also for the lumped mass matrix approaches and therefore, we observe an excellent agreement with the reference solution.  At this point, a few conclusion can be drawn regarding the performance of the proposed SBE. For a moderate number of nodes, ranging from 2 to 9, there is no significant difference in the numerical results regarding the two tested nodal distributions. That is to say, as long as a consistent mass matrix formulation is employed both EQ and GLL points can be recommended and high (even theoretically optimal) ROCs can be achieved. This assessment changes drastically once the mass matrix is diagonalized. In this case, convergence is only achieved for GLL points. Additionally, there was no notable difference between the numerical results obtained by using either the DMM or Guyan reduction techniques. For all analyzed example problems, the results have been virtually identical. A general recommendation is hard to give, since both approaches for treating zero masses have their merits. While the DMM technique is straightforward to implement, it needs an additional parameter to adjust, and is only an approximate method. On the other hand, the Guyan reduction method is exact for the investigated SBEs, but suffers from fill-in of the system matrices, which should not be underestimated for large-scale systems, despite effectively halving the system size.

NUMERICAL EXAMPLES
After having comprehensively discussed the numerical properties of the proposed SBEs, more complex examples are used in the present section to demonstrate the versatility and overall more than satisfactory performance of this element type.
To this end, modal and wave propagation analyses are conducted using tapered beams, where both the material properties, i.e., Young's modulus and the mass density , and the cross-sectional dimensions, i.e., the cross-sectional area and the second moment of area , can be functions of the axial coordinate . The example problems include the free vibration (modal) analysis of a cantilever beam, a cantilever beam subjected to an impact load at the tip, and a damaged beam, where the scattering of the propagating wave packets at the sudden change in cross-section is investigated.

Modal analysis of a tapered beam
In the first numerical example, the eigenfrequencies of a tapered cantilever beam are computed. To this end, we assume a linear variation of the flexural rigidity ( ) as well as a linear variation of the mass per unit length ( ) = ( ) where the parameters 1 and 2 are chosen as 0.95 and 0.8, respectively. That is to say, the flexural rigidity is at the end of the tapered beam is only 5% of the initial value, while the mass per unit length drops to 20% of its initial value. The subscript □ 0 denotes the initial value of a quantity at = 0 (i.e., the left end of the beam). Based on the Frobenius method, an analytical solution for the eigenfrequencies of the tapered beam was derived in Ref. [112], allowing also for the rotation of the beam. The same beam model was also investigated in Ref. [113] by means of the frequency-domain SEM, which is a semi-analytical technique. The results in Refs. [112,113] 9 are provided in terms of a non-dimensionalized frequencyd efined as 9 Remark: In principle, the two functions to describe the variations in the flexural rigidity ( ), which enters the computation of the elemental stiffness matrix in Equation (24), and in the mass per unit length ( ), which enters the computation of the mass matrix in Equation (23), can be chosen arbitrarily. This, however, means that it is not clear what the cross-sectional dimensions of the beam are or how the material properties are changing along the beam axis. The specific values used for the current example are adopted from literature to verify the proposed element. Since both the flexural rigidity and the mass per unit length vary linearly and both 1 and 2 are independent parameters, a change in Young's modulus and the mass density would need to be enforced to obtain the prescribed variations. If a change in geometry of the cross-section is causing the variation in flexural rigidity and mass per unit length, the parameters 1 and 2 would not be independent. In numerical investigations every variation is possible, but in reality the question of manufacturability would be raised.

TA B L E 3
Eigenfrequencies for the tapered cantilever beam. Numerical results are based on 1 spectral beam element with 6 nodes (12 DOFs) and a consistent (CMM) or lumped (LMM) mass matrix formulation.̃, ,̃, , Considering the numerical implementation of the proposed SBE, only minor changes have to be taken into consideration. For each integration point □ , the global coordinate has to be computed by a simple linear mapping where 1 and N represent the nodal coordinates of the first and last node of the beam element, respectively. Thus, the values of ( ) and ( ) can be updated at each integration point. Additionally, the number of integration points for the full-integration is elevated according to the polynomial degree of the spatial variation of the beam properties. This is obviously only done for the elemental stiffness and consistent mass matrices. Since the nodal quadrature technique is utilized to compute the lumped mass matrix, the number of integration points must not be changed. In this particular example, it is sufficient to employ only one additional integration point to accurately account for the linear variation of the properties of the beam during the computation of the elemental matrices.
The tapered beam, analyzed in this section, has the following properties: length = 0.6 m, initial cross-sectional area 0 = 240 × 10 −6 m 2 , initial second moment of area 0 = 2 × 10 −9 m 4 , Young's modulus = 200 GPa, and mass density = 7840kg∕m 3 . Note that the initial values correspond to a rectangular beam of width = 24 mm and height ℎ = 10 mm. To highlight the efficiency of the proposed SBE in conjunction with a consistent mass matrix formulation, the spatial discretization is limited to 1 element with 6 nodes. The results for both the analytical and numerical analysis are listed in Table 3. Here, we observe that the relative error in the first five eigenfrequencies is well below 1%, which is an excellent result considering the coarseness of the mesh. Employing the nodal quadrature technique and using = 1 × 10 7 in the DMM, the results look quite different 10 . While the first three eigenfrequencies are acceptable the fourth and fifth are quite off, which again shows the effect of under-integrating the mass matrix. To obtain results that are of similar accuracy with a lumped mass matrix, the number of elements has to be increased to 3. The numerically obtained values for this case are: 1 = 213.5915, 2 = 972.1761, 3 = 2428.6144, 4 = 4570.9721, and 5 = 7401.5200. Although more DOFs are required, that is, 32 instead of 12 DOFs, the lumped SBEs do still converge and in the end, provides very accurate results. However, there is a lot of potential for further improving the employed mass lumping technique, which is part of ongoing research activities.

Wave propagation analysis in a tapered beam
In order to study the wave propagation in a tapered cantilever beam, a point force is applied at its tip, where the time history of the amplitude functions is selected to approximate a unit impulse excitation-see Section 4.1 for the required input parameters (material properties and dimensions). In Ref. [114], a solution for a rotating beam is presented, which is stiffened due to the emergence of centrifugal forces,. The results have been obtained employing the FDSEM technique. Such stiffening effects are, however, out of the scope of this contribution and therefore, only the non-rotating case is considered. Unfortunately, we are not able to directly compare the numerical results, as information regarding the material damping and the excitation are missing. Nevertheless, a good qualitative agreement is noted. For the analysis, which is presented in the present section, the damping parameters are chosen according to the approach discussed in Appendix A. In this scenario, we assume a damping ratio of 0.01% for the first two modes, which have been already analyzed in Section 4.1. Hence, the values of the damping parameters are R = 0.0350235 and 10 In Section 3.3, we saw that both methods for treating zero masses in the lumped mass matrix formulation yield the same results in modal analyses and therefore, only the results obtained utilizing the DMM technique are discussed at this point. R = 1.68666 × 10 −7 . Furthermore, the unit impulse load is approximated aŝ where 0 denotes the amplitude of the excitation force. The three time parameters, 1 , 2 , and ex represent the times at which the excitation starts, ends, and its duration. Therefore, only two of those time parameters are independent and 2 is defined as For the following transient simulations, the amplitude of the force is 0 = 1 N, the excitation starts at 1 = 1 × 10 −4 s and ends after ex = 1∕22500 s. The resulting forcing function is plotted in Figure 27. For the spatial discretization, we employ only 10 SBEs with 8 nodes, resulting in 142 DOFs. In terms of the temporal discretization, the classical Newmark method (also known as trapezoidal rule or constant average acceleration method) is employed, where a time step size of Δ = 1 × 10 −8 s is selected.
The results of the wave propagation problem in a tapered beam are depicted in Figure 28, where the deflection at the tip of the beam and its velocity are plotted. We observe that all three implementations agree quite well in the recorded displacement history, while distinct differences are noted in the velocity plots. While the solution obtained with adding discrete (rotational) masses, which is referred to as DMM, is in excellent agreement with the results from the consistent mass matrix case, a different velocity response is achieved when employing the Guyan reduction. As discussed in Section 2.3.2, when physical damping is present, the Guyan reduction technique is only an approximate method, even for the proposed SBEs, and therefore, significant errors are expected for high-frequency applications. Thus, even when increasing the number of SBEs and decreasing the time step size for the Newmark method no convergence to the reference results is noted. This constitutes a well-known and severe shortcoming of the Guyan reduction technique when applied transient problems, which is why more advanced model order reduction methods have been developed for transient analyses.
It is very important to note that the two methods for treating zero masses in the lumped mass matrix have only been used to further illustrate their performances. Due to the regularizing effect of the damping matrix on the dynamic stiffness matrix̂, as illustrated by Equation (31), it is actually not required to treat zero masses in the presence of physical damping ( R ≠ 0 and R ≠ 0) in the system. The obtained numerical results are identical to those reported for the DMM approach and therefore, not separately provided for this example.

Damage detection in a tapered beam
In the current section, the wave propagation in a damaged tapered beam, depicted in Figure 29, is studied. To this end, a symmetric notch is introduced by reducing the height of the beam. This notch will cause wave scattering and thus, the beam's response is changed with respect to the healthy state. This problem is of interest for SHM applications, in which guided waves are excited by piezoelectric transducers bonded to the surface of the structure. In the present example, the piezoelectric actuator is modeled by the point-force model [115]. The main idea of this simplified approach is to replace the transducer by two shear forces ( p ( )), which are parallel to the beam axis. The loads are applied in opposite directions at the boundaries of the actuator. Since the nodes of the SBEs are located on/along the beam axis and the point forces need to be directly applied at the nodes, we are required to introduce an additional moment ( p ( ) = p ( )ℎ∕2) to account for the change in the point of application of the concentrated forces. Hence, the relocated point force is responsible for exciting longitudinal waves (not considered in this application) and the moment drives the generation of transverses (flexural) waves. The geometry of the beam is sketched in Figure 29, where the overall length is denoted as and is set to 1.5 m. The cross-section is rectangular with a constant height ℎ of 25 mm and a changing width which follows a quadratic function .
Consequently, the cross-sectional area and the second moment of area also change according to the specified quadratic function for width of the beam. Therefore, two additional integration points are recommended for the computation of the elemental stiffness and consistent mass matrices. The parameters 1 and 2 denote the width at the left and right ends of the beam and are chosen as 100 and 25 mm, respectively. The piezoelectric actuator has a length of 10 mm, while the notch is ℎ∕4 deep (measured from both bottom and top surface, i.e., the thickness of the beam is reduced to ℎ∕2 at the notch) and 10 mm long. Accordingly, the points A, B, C, D, E, and F are located at (0,0), (600,0), (610,0), (1110,0), (1120,0), and (1500,0), respectively. From the known locations of the specific points, the lengths , introduced in Figure 29, can be easily determined and thus, are not explicitly provided at this point. Note that all coordinates are given in the unit [mm], which is omitted for the sake of clarity of notation. The location for F I G U R E 2 9 Sketch of the tapered beam with a notch. recording the response of the beam is referred to as point S, which is located at (900,0). Considering the excitation of the structure, the moments are applied at points B and C with the following amplitude function sin (2 c ).
This function is known as a Hann window modulated sinusoidal tone burst signal, where c is the central frequency and H denotes the duration of the signal, which depends on the number of cycles c and c as H = c ∕ c . In the current example, the structure is excited at 400 kHz with a 4-cycle signal and the beam is made of aluminum, using a Young's modulus of 70 GPa and a mass density of 2700 kg∕m 3 . Note that no Dirichlet boundary conditions have been defined, which is possible due to the regularizing effect of the mass matrix in transient analyses.
The spatial discretization of the structure features 150 SBEs with 5-nodes each resulting in a total of 1202 DOFs. Thus, the element size is 10 mm and hence, the notch is modeled by one element with a reduced height of ℎ∕2. Considering the temporal discretization, Newmark's method (constant average acceleration) is employed with a time step Δ of 1 × 10 −8 s. The results are computed for both the intact and damaged beams (without physical damping) and the recorded wave signals are discussed in the remainder of this section.
In Figure 30, the dynamic response of the tapered beam at the sensor location is shown for both damaged and intact structures. By comparing the results, we clearly observe the influence of the damage, that is, at roughly 3.3 × 10 −5 s an additional wave packet is noted in the damage structure, which is caused by a part of the energy of the traveling wave being reflected at the change in cross-section 11 . After 6.0 × 10 −5 s, the signals obtained from both beams are very different due 11 Remark: According to Ref. [116], the phase and group velocities in an Euler-Bernoulli beam can be determined by p = and g = 2 p . Here, the parameter is related to the material properties and cross-sectional dimensions of the beam as = √ ∕ and the wave number , which can be understood as the "spatial" frequency, is defined as = √ ∕ . Consequently, the phase velocity for an Euler-Bernoulli beam is computed as p = 4 √ ( ∕ ) 2 .
(a) Notch-type damage (b) Intact F I G U R E 3 0 Displacement response recorded at the sensor location for the tapered beam with and without a notch-type damage.
to multiple reflections and scattering of the wave at the notch and at the boundary of the beam. Considering the different numerical methods, we note that a good agreement between the consistent mass matrix approach and both approaches to treat zero masses in the lumped mass matrix is observed. Thus, we can again conclude that for a sufficiently fine spatial and temporal discretization (without the presence of physical damping) both the DMM and Guyan reduction techniques yield very accurate results. Finally, in Figure 31, the wave field throughout the entire (damaged) beam is depicted for different time instances to get a better idea of the behavior of the traveling wave packets. For the sake of clarity, the actuator is visualized as a gray rectangle, while the location of the notch is indicated by a blue rectangle. Additionally, the sensor location is shown as a red filled circle. In Figure 31a, the wave field is plotted at 1 = 1 × 10 −5 s. At this point, the excitation is finished and therefore, we only observe the initial wave packets traveling in both directions. A reflected wave packet originating from the notch is observed in Figure 31b at 2 = 4 × 10 −5 s, while the initially excited wave already reaches the left boundary of the beam and interacts with it. The displacement field at 3 = 6 × 10 −5 s is plotted in Figure 31c, where we notice wave reflections at the right boundary of the tapered beam. Due to the multiple interactions of the wave packets with the notchtype damage and the boundaries of the beam, the complexity of the measures signal is significantly increased and we observe the distinct differences between the intact and damaged structures as previously indicated in Figure 30.  By means of this more complex example, we also want to illustrate the advantages of using high-order spectral beam elements in comparison to standard 2-noded beam elements. To this end, the simulation is repeated for spatial discretizations employing 150 beam elements featuring 2, 3, 4, 5, and 6 nodes (p-refinement). Considering an h-refinement strategy, 2noded beam elements are deployed, where the structure is divided into 150, 300, 600, 1200, and 2400 elements, respectively. In order to exclude effects stemming from the time integration scheme, the time step size is reduced to Δ = 1 × 10 −10 s, while still applying Newmark's constant average acceleration method (trapezoidal rule). In Figure 32, the obtained convergence curves are plotted over the number of degrees of freedom (NDOF) as wells as over the computational time. From the results depicted in Figure 32, we can conclude that significant savings in terms of NDOF are possible when utilizing SBEs, which directly translate to savings in the computational time. It can be easily observed that to obtain results with similar accuracy, the use of high-order elements requires much less DOFs as well as a shorter computational time. Note that the error has been computed in the 2 -norm where the reference values are obtained numerical by means of an overkill solution using 300 SBEs with 8 nodes each, resulting in a high-fidelity analysis.

CONCLUSIONS
In the article at hand, a thorough numerical analysis of the properties of the previously introduced SBEs [74] is conducted. Overall, the findings are very encouraging but ongoing research on suitable mass lumping methods for 1 -continuous elements is still of paramount importance. In the following, the main findings are summarized and important conclusions are drawn. It has been found that: 4. Both methods for treating zero masses after lumping, i.e., DMM and Guyan reduction, achieve almost identical results if no physical damping is introduced. 5. In cases where physical damping is of importance, no method to treat zero masses is required due to the regularizing effect of the damping matrix , provided that both Rayleigh parameters R and R are different from zero. On the other hand, if the physical damping is only mass-proportional, i.e., R ≠ 0 and R = 0, the above conclusions are still valid. 6. In terms of the selected nodal distribution to derive the shape functions, no significant differences in the results are observed for a moderate number of nodes and therefore, both EQ and GLL points can be employed without reservation. 7. To achieve converged results for a lumped mass matrix formulation, GLL points have to be used, due to the better accuracy in the integration of the mass matrix when using GLL quadrature rules compared to Newton-Cotes formulas, which are applicable for elements based on EQ points.
In a nutshell, the investigated SBEs offers highly accurate results in both static and dynamic analyses, especially in combination with a consistent mass matrix formulation. Due to the formulation of the element's shape functions, mass lumping is possible and a viable option for explicit dynamics. However, ongoing research must be directed at developing novel mass lumping methods to effectively diagonalize the mass matrix of beam elements, guaranteeing high-order convergence rates.
Finally, we want to conclude with a remark on the extension of the presented SBEs to plate or shell problems. In this context, the reader is reminded of the classical Bogner-Fox-Schmit (BFS) element [117,118], which is, at least in our opinion, the simplest formulation (also compared to IGA) to attain full 1 -continuity along the element edges. BFS elements are straightforwardly derived by a simple tensor product extension of our SBE shape functions, meaning that merely the proposed one-dimensional Hermitian shape functions are utilized to obtain the two-dimensional shape functions (in the reference domain), which are required for plate and shell problems. To counteract an inherent limitation related to BFS elements, which is the requirement of employing a regular/structured mesh 12 , a combination with fictitious or embedded domain methods is advised. This is fully in line with recent research activities on reviving BFS elements as proposed in Refs. [11,12], where static problems or linear buckling analyses have been tackled. In this regard, it is expected that the numerical properties determined in this paper are also exhibited by high-order cut BFS elements based on a tensor product formulation of the one-dimensional shape functions.

A C K N O W L E D G M E N T S
The authors have nothing to report.
Open access funding enabled and organized by Projekt DEAL.

APPENDIX A: RAYLEIGH DAMPING
In all dynamical applications, damping is a highly problematic topic, which still needs dedicated research efforts. In our implementation, we simply rely on viscous damping and utilize Rayleigh's hypothesis, which states that the damping matrix is a linear combination of the stiffness and mass matrices = R + R , (A.1) where R and R denote the mass-and stiffness-proportional damping coefficients. These two parameters are commonly determined by prescribing certain damping ratios for two distinct modes of vibration. This, however, can hardly capture all effects of a complex structure, since usually more than just two modes contribute in a significant manner to the dynamic behavior of a structure. Although, this constitutes a severe shortcoming of Rayleigh's hypothesis, it is nonetheless the most widely-used approach to account for material damping. Let us recall the conventional approach to identifying the damping coefficients R and R . To this end, we select two mode shapes that significantly contribute to the dynamic behavior and prescribe their damping ratios , which must be determined by experiments before conducting simulations. Together with the natural frequencies of the selected mode shapes , we can calculate the damping coefficients by solving the following system of equations Note that the modes of vibration have to be chosen very carefully as otherwise the transient response might be significantly distorted and thus, the numerical response is not representative of the physics of the system. In Figure A.1, the variation of the damping ratio with respect to the circular frequency is depicted for three different cases. The solid black line ( : ) represents the complete Rayleigh damping, while the (thin) solid grey lines illustrate the damping behavior of the system for mass-( m : ) or stiffness-proportional damping ( k : ), respectively. We clearly observe, that the damping ratio is exactly prescribed for the two selected frequencies, while all frequencies within this interval exhibit a smaller damping ratio, which is still reasonably close to the prescribed values. Outside of the specified frequency range, the damping ratio increases monotonously without bounds.
In Ref. [119], a novel methodology was introduced to improve the computation of the damping coefficients R and R , with the aim of extending the simple Rayleigh damping technique to more complex structures. Based on the use of mode participation factors (masses), the dominant vibration modes are selected. Due to the fact that an arbitrary number of modes can be selected, a least squares fit is proposed to determine the two damping coefficients. This strategy results in an improved agreement between experiments and numerical simulations. The authors additionally provide a good review and discuss other possible strategies to tackle the shortcomings of Rayleigh's hypothesis. In Ref. [120], an extended Rayleigh damping model is proposed to increase the nearly constant damping range mentioned above. To this end, the stiffness-proportional damping part is exchanged with a causal damping model, involving more constants that need to be determined. On the other hand, in Ref. [121] a methodology is presented to adjust the numerical damping to the physical one. Based on the natural frequencies and modal damping parameters of a structure, which are determined experimentally, a generalized proportional damping technique is developed, where the damping matrix is expressed as a function of the product −1 . In principle, arbitrary variations of the damping ratio with the frequency can be taken into account, but matrix functions are rather costly to compute for large-scale systems and the derived damping matrix is fully-populated. Despite the mentioned research efforts, there is certainly still a lot of research to be done to accurately model the physical damping of complex structures. So far, satisfactory solutions exist only for particular problems, but their extension to a general framework is hardly ever possible. Therefore, we rely on Rayleigh's hypothesis for the numerical analysis of the proposed SBE, as it is unfortunately out of the scope of the current article to describe the damping behavior of complex structures in a comprehensive way.