Isogeometric Analysis of Structural Vibrations — Source link

This paper begins with personal recollections of John H. Argyris. The geometrical spirit embodied in Argyris’s work is revived in the sequel in applying the newly developed concept of isogeometric analysis to structural vibration problems. After reviewing some funda-mentals of isogeometric analysis, application is made to several structural models, including rods, thin beams, membranes, and thin plates. Rotationless beam and plate models are utilized as well as three-dimensional solid models. The concept of k -reﬁnement is explored and shown to produce more accurate and robust results than corresponding ﬁnite elements. Through the use of nonlinear parameterization, ‘‘optical’’ branches of frequency spectra are eliminated for k -reﬁned meshes. Optical branches have been identiﬁed as contributors to Gibbs phenomena in wave propagation problems and the cause of rapid degradation of higher modes in p -method ﬁnite elements. A geometrically exact model of the NASA Aluminum Testbed Cylinder is constructed and frequencies and mode shapes are computed and shown to compare favorably with experimental results.


Thomas J.R. Hughes's personal recollections of John H. Argyris
John Argyris was my initial inspiration in science. I first heard of the finite element method in 1967 while I was working at General Dynamics, Electric Boat Division in Groton, Connecticut. Scientific journals were made available to staff and one, in particular, caught my attention, the Aeronautical Journal of the Royal Society which contained short articles in each issue by John Argyris reporting the latest developments in the finite element method. I seem to recall that the terminologies ''matrix displacement method'' and ''matrix force method'' were still given preference, but eventually they faded away along with other synonyms for finite elements. I remember the meticulously prepared, elaborate drawings in those papers, a hallmark of John's scientific work. John was at the peak of his scientific powers, his writing was eloquent, and he was clearly having fun developing and naming new elements and element families, such as ''Lumina,'' ''Hermes,'' ''Sheba,'' etc. I was trying to understand what was going on in this revolutionary new field. A frequent reference in these works was Argyris and Kelsey's Energy Theorems and Structural Analysis [3], a republication of a series of articles which originally appeared in the journal Aircraft Engineering in the mid 1950s [2]. I ordered this book and I was again fascinated by the carefully done illustrations of what seemed to me to be very complex structural geometries, certainly more complex than the ones appearing in structural analysis texts of that time. One of the articles in this series contained the famous diptych

Introduction
In Hughes et al. [27] we introduced the concept of isogeometric analysis, which may be viewed as a logical extension of finite element analysis. The objectives of the isogeometric approach were to develop an analysis framework based on functions employed in computer aided design (CAD) systems, capable of representing many engineering geometries exactly; to employ one, and only one, geometric description for all meshes and all orders of approximation; and to vastly simplify mesh refinement procedures. As a primary tool in the establishment of this new framework for analysis, we selected NURBS (non-uniform rational B-Splines; see, e.g., [38,35]). We found NURBS to possess many interesting properties in analysis and excellent results were attained for problems of linear solid and structural mechanics and linear shells modeled as three-dimensional solids, and particularly intriguing results were obtained for the linear advection-diffusion equation. Indeed, it was concluded that isogeometric analysis provided a viable alternative to existing finite element analysis procedures and possessed a number of advantages that might be exploited in various situations. However, isogeometric analysis is in its infancy and much basic work remains to be done to bring these ideas to fruition.
A fundamental tenet of isogeometric analysis is to represent geometry as accurately as possible. It was argued in [27] that the faceted nature of finite element geometries could lead to significant errors and difficulties. This is schematically conveyed in Fig. 1. In order to generate meshes, geometrical simplifications are introduced in finite element analysis. For example, features such as small holes and fillets are often removed. Stress concentrations produced by holes are then missing, and artificial, non-physical, stress concentrations are induced by the removal of fillets. The stresses at sharp, reentrant corners will be infinite, which makes adaptive mesh refinement strategies meaningless. If the refinement is performed to capture geometrical features in the limit, then tight, automated communication with the geometry definition, typically a CAD file, must exist for the mesh generator and solver. It is rarely the case that this ideal situation is attained in industrial settings, which seems to be the reason that automatic, adaptive, refinement procedures have had little industrial penetration despite enormous academic research activity. In isogeometric analysis, the first mesh is designed to represent the exact geometry, and subsequent refinements are obtained without further communication with the CAD representation. This idea is dramatized in Fig. 2 in which the question ''What is a circle?'' is asked rhetorically. In finite element analysis, a circle is an ideal achieved in the limit of mesh refinement (i.e., h-refinement) but never achieved in reality, whereas a circle is achieved exactly for the coarsest mesh in isogeometric analysis, and this exact geometry, and its parameterization, are maintained for all mesh refinements. It is interesting to note that, in the limit, the isogeometric model converges to a polynomial representation on each element, but not for any finite mesh. This is the obverse of finite element analysis in which polynomial approximations exist on all meshes, and the circle is the idealized limit.
In this paper we initiate the study of the isogeometric analysis methodology in structural vibration analysis. In Section 3 we briefly review the basic concepts of isogeometric analysis. The interested reader is also urged to consult our recent work [27] which presents a more comprehensive introduction. We emphasize the concept of k-refinement, a higher-order procedure employing smooth basis functions, which is used repeatedly in the vibration calculations later on. In Sections 4-7 we investigate isogeometric approaches to some simple model problems of structural vibration, including the longitudinal vibration of a rod (and, equivalently, the transverse vibration of a string, or shear beam), the transverse vibration of a thin beam governed by Bernoulli-Euler theory, the transverse vibration of membranes, the transverse vibration of thin plates governed by Poisson-Kirchhoff theory, and the transverse vibration of a thin plate modeled as a three-dimensional elastic solid. In the cases of Bernoulli-Euler beam theory and Poisson-Kirchhoff plate theory, we have employed rotationless formulations, an important theme of contemporary research in structural mechanics (see, e.g., [14][15][16]19,[31][32][33][34]). 2 In the one-dimensional cases we performed numerical analyses of discrete frequency spectra. We were also able to theoretically derive the continuous, limiting spectra and we determined that these spectra are invariant if normalized by the total number of degrees-of-freedom in the model. In other words, one is able to determine a priori the error in frequency for a particular mode from a single function, no matter how many degrees-of-freedom are present in the model. These elementary results are very useful in determining the vibration characteristics of isogeometric models and provide a basis for comparison with standard finite element discretizations. It is well known that, in the case of higher-order finite elements, ''optical'' branches are present in the spectra (see [8]) and that these are responsible for the large errors in the high-frequency part of the spectrum (see [26]) and contribute to the oscillations (i.e., ''Gibbs phenomena'') that appear about discontinuities in wave propagation problems. The accurate branch, the so-called ''acoustic'' branch (see [8]), corresponds to the low-frequency part of the spectrum. In finite element analysis, both acoustic and optical branches are continuous, and the optical branches vitiate a significant portion of the spectrum. In isogeometric analysis, when a linear parameterization of the geometrical mapping from the patch to its image in physical space is employed, only a finite number of frequencies constitute the optical branch. The number of modes comprising the optical branch is constant once the order of approximation is set, independent of the number of elements, but increases with order. In this case, almost the entire spectrum corresponds to the acoustic branch. A linear parameterization of the mapping requires a non-uniform distribution of control points. Our previous work [27] describes the algorithm which locates control points to attain a linear parameterization. Spacing control points uniformly produces a nonlinear parameterization of the mapping. In this case, remarkably, the optical branch is entirely eliminated! The convergence rates of higher-order finite elements and isogeometric elements constructed by k-refinement are the same for the same order basis, but the overall accuracy of the spectrum is much greater for isogeometric elements. These observations corroborate the speculation that the k-method would be a more accurate and economical procedure than p-method finite elements in vibration analysis of structural members. Studies of membranes and thin plates provide additional corroboration. We also present some initial studies of mass lumping within the isogeometric approach. The ''row sum'' technique is employed (see [26]). Due to the pointwise non-negativity of B-Spline and NURBS bases, the row sum technique is guaranteed to produce positive lumped masses but only second-order accurate frequencies are obtained, independently of the order of basis functions employed. This is unsatisfactory but we conjecture that, by appropriately locating knots and control points, higher-order-accurate lumping procedures may exist. This is a topic requiring further research.
In Section 8, we apply the isogeometric approach to the NASA Aluminum Testbed Cylinder (ATC) which has been extensively studied experimentally to determine its vibration characteristics. Our isogeometric model is an exact threedimensional version of the ''as-drawn'' geometry. All fine-scale features of the geometry, such as fillets, are precisely accounted for. Comparisons are made between the experimental data and the numerical results.
In Section 9 we draw conclusions. Appendix A presents analytical and numerical results concerning the order of accuracy of consistent and lumped mass schemes.

A brief summary of NURBS-based isogeometric analysis
Non-uniform rational B-Splines (NURBS) are a standard tool for describing and modeling curves and surfaces in computer aided design and computer graphics (see [35,38]). The aim of this section is to introduce them briefly and to present an overview of isogeometric analysis, for which an extensive account has been given in Hughes et al. [27]. ''What is a circle?'' In finite element analysis it is an idealization attained in the limit of mesh refinement but never for any finite mesh. In isogeometric analysis, the same exact geometry and parameterization are maintained for all meshes.

B-Splines
B-Splines are piecewise polynomial curves composed of linear combinations of B-Spline basis functions. The coefficients are points in space, referred to as control points.

Knot vectors
A knot vector, N, is a set of non-decreasing real numbers representing coordinates in the parametric space of the curve: where p is the order of the B-Spline and n is the number of basis functions (and control points) necessary to describe it. The interval [n 1 , n n+p+1 ] is called a patch. A knot vector is said to be uniform if its knots are uniformly spaced and non-uniform otherwise. Moreover, a knot vector is said to be open if its first and last knots are repeated p + 1 times. In what follows, we always employ open knot vectors. Basis functions formed from open knot vectors are interpolatory at the ends of the parametric interval [n 1 , n n+p+1 ] but are not, in general, interpolatory at interior knots.

Basis functions
Given a knot vector, N, B-Spline basis functions are defined recursively starting with p = 0 (piecewise constants): For p >1: In Fig. 3 we present an example consisting of n = 9 cubic basis functions generated from the open knot vector N = {0, 0, 0, 0, 1/6, 1/3, 1/2, 2/3, 5/6, 1, 1, 1, 1}. An important property of B-Spline basis functions is that they are C pÀ1 -continuous, if internal knots are not repeated. If a knot has multiplicity k, the basis is C pÀk -continuous at that knot. In particular, when a knot has multiplicity p, the basis is C 0 and interpolatory at that location.
Other noteworthy properties are • B-Spline basis functions formed from an open knot vector constitute a partition of unity, that is, P n i¼1 N i;p ðnÞ¼1 8n. • The support of each N i,p is compact and contained in the interval [n i , n i+p+1 ].
• B-Spline basis functions are non-negative: N i,p (n) P 0 "n.

B-Spline curves
Given the order of a desired B-Spline and an appropriately defined knot vector, it is possible to construct n basis functions. In turn, given a set of n control points in R d , the piecewise polynomial B-Spline curve C(n) of order p is obtained by taking a linear combination of basis functions and control points: where B i is the ith control point. The piecewise linear interpolation of the control points defines the control net.InFig. 4 we present a cubic two-dimensional B-Spline curve generated from the basis functions shown in Fig. 3 along with its control net. A B-Spline curve has continuous derivatives of order p À 1, which can be decreased by k if a knot or a control point has multiplicity k + 1. An important property of B-Spline curves is affine covariance, which states that an affine transformation of the curve is obtained by applying the transformation to its control points.

B-Spline surfaces
By means of tensor products, B-Spline surfaces can be constructed starting from knot vectors N ={n 1 , ..., n n+p+1 } and H ={g 1 , ..., g m+q+1 }, and an n · m net of control points B i,j . One-dimensional basis functions N i,p and M j,q (with i =1,..., n and j =1,..., m) of order p and q, respectively, are defined from the knot vectors, and the B-Spline surface is constructed as The patch is the domain [n 1 ,n n+p+1 ] · [g 1 ,g m+q+1 ].

Non-uniform rational B-Splines
A rational B-Spline in R d is the projection onto d-dimensional physical space of a polynomial B-Spline defined in (d + 1)-dimensional homogeneous coordinate space. For a complete discussion of these space projections, see [20] and references therein. In this way, a great variety of geometrical entities can be constructed and, in particular, all conic sections can be obtained exactly. The projective transformation of a B-Spline curve yields a rational polynomial curve. Note that when we refer to the ''order'' of a NURBS curve, we mean the order of the polynomial curve from which the rational curve was generated.
To obtain a NURBS curve in R d , we start from a set B w i (i =1,..., n) of control points (''projective points'') for a B-Spline curve in R dþ1 with knot vector N. Then the control points for the NURBS curve are where (B i ) j is the jth component of the vector B i and w i ¼ðB w i Þ dþ1 is referred to as the ith weight. The NURBS basis functions of order p are then defined as and their first and second derivatives are given by and ð11Þ The NURBS curve is defined by Rational surfaces and solids are defined in an analogous way in terms of the basis functions (resp.) and R p;q;r i;j;k ðn; g; fÞ¼ In the following, we summarize noteworthy properties of NURBS: • NURBS basis functions formed from an open knot vector constitute a partition of unity: P n i¼1 R p i ðnÞ¼1 8n. • The continuity and supports of NURBS basis functions are the same as for B-Splines.
• NURBS possess the property of affine covariance.
• If all weights are equal, NURBS become B-Splines.
• NURBS surfaces and solids are the projective transformations of tensor product piecewise polynomial entities.

Isogeometric analysis
A summary of the main features of isogeometric analysis follows: • A mesh for a NURBS patch is defined by the product of knot vectors. For example, in three-dimensions, a mesh is given by N · H · Z. • Knot spans subdivide the domain into ''elements.'' • The support of each basis function consists of a small number of elements.
• The control points associated with the basis functions define the geometry.
• The isoparametric concept is invoked, that is, the unknown variables are represented in terms of the basis functions which define the geometry. The coefficients of the basis functions are the degrees-of-freedom, or control variables. • Three different mesh refinement strategies are possible: analogues of classical h-refinement (by knot insertion) and prefinement (by order elevation of the basis functions), and a new possibility referred to as k-refinement, which increases smoothness in addition to order. • The element arrays constructed from isoparametric NURBS can be assembled into global arrays in the same way as finite elements (see [26,Chapter 2]). Compatibility of NURBS patches is attained by employing the same NURBS edge and surface representations on both sides of patch interfaces. This gives rise to a standard continuous Galerkin method and mesh refinement necessarily propagates from patch to patch. There is also the possibility of employing discontinuous Galerkin methods along patch boundaries.
• Dirichlet boundary conditions are applied to the control variables. If the Dirichlet conditions are homogeneous, this results in exact pointwise satisfaction. If they are inhomogeneous, the boundary values must be approximated by functions lying within the NURBS space, and this results in ''strong'' but approximate satisfaction of the boundary conditions, as in finite elements. Another option is to impose Dirichlet conditions ''weakly'' (we will discuss this later on). Neumann boundary conditions are satisfied naturally as in standard finite element formulations (see [26,Chapters 1 and 2]).
In structural analysis, NURBS elements represent all rigid body motions and constant strain states exactly (see [26]). Consequently, structures assembled from compatible NURBS elements pass standard ''patch tests'' (see [26,Chapters 3 and 4], for a description of patch tests).

k-Refinement
Isogeometric analysis is fundamentally a higher-order approach. While it is true that the first two NURBS bases consist of constants and linears, identical in every way to standard finite elements, it takes at least quadratic-level NURBS to exactly represent conic sections. Refinement procedures are also fundamental components of NURBS technology. There are analogues of finite element h-a n dp-refinement procedures, and a new, potentially more efficient, higher-order procedure, k-refinement (see [27]). In p-refinement, C 0 -continuity is maintained across knots (i.e., ''element'' boundaries). In krefinement, continuity of order C pÀ1 is attained across knots, at least within patches. The additional smoothness in k-refinement seems intuitively appealing for situations in which exact solutions are dominantly very smooth, such as free vibrations of structures and bifurcation buckling of thin beams, plates and shells. k-refinement also offers a very concise parameter- ization of smooth functions. The potential efficiency gains of k-refinement are suggested by the following calculations comparing p-and k-refinement. First, consider a one-dimensional mesh with n basis functions of order p. Note that the number of basis functions is equal to the number of control variables, and is also equal to the number of control points. After r refinements (i.e., order elevations), the number of basis functions, each of order p + r,is(r +1)n À rp for p-refinement and n + r for k-refinement. The growth in the number of control variables is depicted graphically in Fig. 5. Next, consider a ddimensional mesh with n d basis functions. After r refinements, assuming r to be large, the number of basis functions asymptotically approaches n d (r d + dr dÀ1 ) for p-refinement and n d (1 + drn À1 ) for k-refinement. The difference is seen to be very significant. Graphical comparisons for two and three dimensions are presented in Figs. 6 and 7, respectively. Keep in mind that the mesh, defined by the knot locations, is fixed and is the same for p-and k-refinement. See Hughes et al. [27] for further details.

Natural frequencies and modes
Given a linear multi-degree-of-freedom structural system, the undamped, unforced equations of motion, which govern free vibrations, are where M and K are, respectively, the consistent mass and the stiffness matrices, u = u(t) is the displacement vector and € u ¼ d 2 u dt 2 is the acceleration vector. Due to the fact that B-Spline and NURBS bases are pointwise non-negative, every entry of the consistent mass matrix is non-negative. The nth normal mode, / n , can be obtained by separation of variables: where x n is the nth natural frequency. Combining (15) and (16) leads to the generalized eigenproblem: The normal modes are defined up to a multiplicative constant. Different ways of normalization have been proposed. One of the most widely used is

One-dimensional problems
In the following examples, due to the simplicity of the geometry, all of the weights are equal to 1 (i.e., NURBS basis functions become B-Splines).

Longitudinal vibrations of an elastic rod
We study the problem of the structural vibrations of an elastic fixed-fixed rod of unit length, whose natural frequencies and modes, assuming unit material parameters, are governed by and for which the exact solution in terms of natural frequencies is After writing the weak formulation and performing the discretization, a problem of the form of (17) is obtained.

Numerical experiments
As a first numerical experiment, the generalized eigenproblem (17) is solved with both finite elements and isogeometric analysis using quadratic basis functions. The resulting natural frequencies, x h n , are presented in Fig. 8, normalized with respect to the exact solution (20), and plotted versus the mode number, n, normalized by the total number of degreesof-freedom, N. To produce the spectra of Fig. 8, we used N = 999 but the results are in fact independent of N. Fig. 8 illustrates the superior behavior of NURBS basis functions compared with finite elements. In this case, the finite element results depict an acoustical branch for n/N < 0.5 and an optical branch for n/N > 0.5 (see [8]).
We then performed the same eigenvalue analysis using higher-order NURBS basis functions. The resulting spectra are presented in Fig. 9; the analyses were carried out using N = 1000 degrees-of-freedom.
Increasing the order, p, of the basis functions, the results show higher accuracy, namely, 2p (see Appendix A for the computation of the order of accuracy using quadratic and cubic NURBS). Figs. 11-13 confirm that the order of convergence for frequencies computed using NURBS is O(h 2p ), as with polynomial-based finite elements. Increasing p also results in the appearance of strange frequencies at the very end of the spectrum (see Fig. 9), referred to in the following as ''outlier frequencies'' (in analogy with outlier values in statistics, see, e.g., [30]), whose number and magnitude increase with p.I n Fig. 10, this behavior is highlighted by plotting the last computed frequencies for p =2,..., 10. To understand the outliers, we first remark that the finite element spectrum for quadratic elements consists of acoustic (low-mode) and optical (highmode) branches, in the sense of [8]. Both these branches are continuous as may be seen from Fig. 8. There are only two distinct equations in the discrete system, corresponding to element middle and end nodes, and this gives rise to the two branches. In the case of NURBS, all but a finite number of equation are the same. The ones associated with the open knot vectors are different, and are responsible for the outliers. The outliers constitute a discrete optical branch. The typical equation of the interior knots gives rise to the continuous acoustic branch, as will be analytically verified in the next section. In finite element analysis, the frequencies associated with the optical branch are regarded as inaccurate and, obviously, the same is true for NURBS. In many applications, these frequencies are harmless. They can be ignored in vibration analysis and their participation in transient response can be suppressed through the use of dissipative, implicit time integration algorithms (see, e.g., [13,29,[24][25][26]). However, they would be detrimental in explicit transient analysis because the frequencies of   the highest modes are grossly overestimated and stability would necessitate an unacceptably small time step, but it will be shown in the next section how to completely eliminate the outliers by a reparameterization of the isogeometric mapping.

Analytical determination of the discrete spectrum
Following the derivations of Hughes [26,Chapter 9], it is possible to analytically compute the discrete spectra previously determined numerically. The starting point is the mass and stiffness matrices for a typical interior element (note that for interior elements, the basis functions are all identical). For quadratic NURBS, the mass and stiffness of a typical interior element are where h =1/n el =1/(N À p), n el is the number of elements, N is the number of control points, and p = 2 is the order of the basis functions. The equation of motion for the typical interior control point, A,i s  which can be compactly written as where a and b are operators defined as follows: Separating variables and substituting this expression into (23), after adding and subtracting ðx h hÞ 2 20 au A , we obtain: The satisfaction of (26) is achieved by selecting / A and q such that: Assuming a solution for (27) of the form (for fixed-fixed boundary conditions): Eq. (27) can be rewritten as Substituting expressions (24) for a and b, and using the trigonometric identity sinða AE bÞ¼sinðaÞ cosðbÞAEsinðbÞ cosðaÞ, yields: ðx h hÞ 2 20 ð16 þ 13 cosðxhÞþcos 2 ðxhÞÞÀð2 À cosðxhÞÀcos 2 ðxhÞÞ ¼ 0; which can be solved for x h x , giving Eq. (32) is the analytical expression for the normalized discrete spectrum for our problem, using quadratic NURBS basis functions. Analogous calculations can be performed for higher-order approximations. The expression for cubic NURBS, is In Fig. 14 we present the analytical and numerical spectra for quadratic and cubic NURBS. For the computation of the numerical spectra, 2000 control points were employed. The only differences are the outlier frequencies at the end of the numerical spectrum obtained for cubic NURBS. Remark 2. All the numerical results described up to now have been obtained using control points computed with the procedure proposed by Hughes et al. [27], which leads to linear parameterization (i.e., constant Jacobian determinant).
The results obtained are seen to be very good, except for the outliers, which get progressively worse for higher-order approximations. A way to avoid this behavior is to employ uniformly spaced control points. The difference between a distribution of 21 control points in the case of linear parameterization and of uniformly spaced points, using cubic NURBS, is presented in Fig. 15. This choice corresponds to a nonlinear parameterization (see Figs. 16 and 17 for plots of the parameterization x(n) and its Jacobian J ðnÞ¼ dxðnÞ dn for the cases in Fig. 15). Fig. 18 presents spectra computed using uniformly spaced control points. The outlier frequencies are eliminated and the continuous spectra coincide with the ones computed analytically and presented previously in Fig. 14    Remark 3. In this work, consistent mass is emphasized because it seems more suitable than lumped mass when higherorder approximations are involved. However, some preliminary tests were performed with a ''row sum'' lumped mass (see [26,Chapter 7]). This approach proves satisfactory for some low-order finite elements but it is incapable of maintaining full accuracy in the present context (see Fig. 19). In all cases, accuracy is limited to second order. Analytical and numerical lumped mass results for quadratic and cubic NURBS are presented in Appendix A. Despite these negative results, we do not think the issue of lumped mass and NURBS is closed. There may be ways to develop higher-order accurate lumped mass matrices. Inspiration may be taken from the work of Fried and Malkus [21]. Perhaps nonlinear parameterizations and nonuniform knot distributions, in conjunction with a lumping scheme, are worthwhile directions to explore. This seems to be an interesting problem of applied mathematics with practical significance.

Transverse vibrations of an Bernoulli-Euler beam
The transverse vibrations of a simply supported, unit length Bernoulli-Euler beam are considered (see [26,Chapter 7]). For this case, the natural frequencies and modes, assuming unit material and cross-sectional parameters, are governed by u ;xxxx À x 2 u ¼ 0 for x 20; 1½; uð0Þ¼uð1Þ¼u ;xx ð0Þ¼u ;xx ð1Þ¼0; ð34Þ

Numerical experiments
The numerical experiments and results for the Bernoulli-Euler beam problem are analogous to the ones reported for the rod. A remark about the formulation is in order before presenting the results. The classical beam finite element employed to solve problem (34) is a two-node Hermite cubic element with two degrees-of-freedom per node (transverse displacement and rotation), whereas our isogeometric analysis formulation is rotation-free (see, for example, Engel et al. [19]). Later in this section we will discuss the problem of the imposition of rotation boundary conditions. Discrete spectra obtained using classical cubic finite element and NURBS basis functions are presented in Fig. 20. The NURBS solution is significantly more accurate but two outlier frequencies are present at the end of the spectrum. Fig. 21 presents the discrete spectra obtained using different order NURBS basis functions. The behavior is similar to the case of the rod, including the outlier frequencies (see Fig. 22). Note that quadratic NURBS are admissible in the present context because they are C 1 -continuous on patches. Slope continuity may be weakly enforced across patch boundaries by way of the technique described in Engel et al. [19]. There are no outliers for quadratic NURBS but the accuracy level is rather poor compared with cubics.
Figs. [23][24][25] show that the order of convergence of frequencies using NURBS is optimal, that is O(h 2(pÀ1) ).   The analytical computation of the discrete spectrum, performed previously for the rod problem, can also be done in the present case. Employing cubic NURBS shape functions, for example, gives rise to the following expression: The analytical and numerical discrete spectra for cubic and quartic approximations are compared in Fig. 26. For the computation of the numerical discrete spectra, 2000 control points were used. The only differences are in the outlier frequencies at the end of the numerical discrete spectra. As previously shown for the rod problem, the outliers can be removed by nonlinear parameterization derived from a uniformly spaced distribution of control points. In this way, the discrete spectra of Fig. 27 are obtained, which coincide with the analytically computed ones.

Rotation boundary conditions
The Bernoulli-Euler beam formulation employed is ''rotation-free,'' that is, only displacements are degrees-of-freedom. Rotations (i.e., slopes) can be computed as derivatives of displacement but are not degrees-of-freedom. To illustrate the method utilized to enforce rotation boundary conditions, we consider the following problem of a cantilever beam: The natural frequencies are (see [12]) x n ¼ b 2 n ; with b 1 = 1.8751, b 2 = 4.6941, b 3 = 7.8548, b 4 = 10.996, and b n =(n À 1/2)p for n > 4. Two strategies were employed to solve this problem. One is based on weak boundary condition imposition and the other on Lagrange multipliers. The former is the approach used in Engel et al. [19]. In this case, the bilinear form, from which the stiffness matrix derives, is given by where v h and u h are the discrete weighting and trial solution, respectively, and s is a stabilization parameter. Analogous to what is done in Prudhomme et al. [37] for the Poisson problem, it can be shown that the choice of s needs to be proportional to p 2 /h, where p is the order of the NURBS basis and h is the mesh parameter. With this formulation, the cantilever beam problem (37) was solved, and corresponding discrete spectra for different order NURBS are shown in Fig. 28 (1000 control points and s = p 2 /h were used). Fig. 28 shows the same behavior seen in Fig. 21 for the simply supported case. The other way to enforce rotation boundary conditions is through Lagrange multipliers. In this case the bilinear form is  where k is the Lagrange multiplier and l is its weighting counterpart. The advantage of the Lagrange multiplier approach is that the rotation boundary condition is exactly enforced. Results for the Lagrange multiplier approach are presented in Fig. 29. For all practical purposes, the results of the two approaches are the same (cf., Figs. 28 and 29).

Two-dimensional problems
In this section we present some numerical experiments for two-dimensional counterparts of the rod and Bernoulli-Euler beam problems considered previously, namely, the transverse vibrations of an elastic membrane and transverse vibrations of a Poisson-Kirchhoff plate, respectively.

Transverse vibrations of an elastic membrane
The first problem we consider consists of the study of the transverse vibrations of a square, elastic membrane, whose natural frequencies and modes, assuming unit tension, density and edge length, are governed by  where $ 2 is the Laplace operator. The exact natural frequencies are (see, e.g., [28]): The numerical results are qualitatively similar to the ones obtained in the study of the one-dimensional problems. The normalized discrete spectra obtained employing different order NURBS basis functions and using a linear parameterization over a 90 · 90 control net are presented in Fig. 30. Note that l is the number of modes sorted from the lowest to the highest in frequency, while N is the total number of degrees-of-freedom. Fig. 31 shows the lower half of the frequency spectra to highlight the accuracy of the different approximations. The optical branches seen in Fig. 30 can again be eliminated by a uniformly spaced control net, as shown in Fig. 32.

Transverse vibrations of a Poisson-Kirchhoff plate
We consider the transverse vibrations of a simply supported, square plate governed by Poisson-Kirchhoff plate theory. The natural frequencies and modes, assuming unit flexural stiffness, density and edge length, are governed by the biharmonic problem: r 4 uðx; yÞÀx 2 uðx; yÞ¼0 for ðx; yÞ2X ¼0; 1½Â0; 1½; uðx; yÞj oX ¼ 0; r 2 uðx; yÞj oX ¼ 0; for which the exact solution natural frequencies (see, e.g., [28]) are For this case, as for the Bernoulli-Euler beam, the NURBS formulation results in a rotation-free approach. The boundary conditions on rotations can be imposed in similar fashion to the way described for the beam (see [19] for further details). The numerical results are similar to the ones obtained for the elastic membrane. In Fig. 33, the normalized discrete spectra using a linear parameterization and a 90 · 90 control net are presented. Fig. 34 shows a detail of the lower-frequency part. The y-axis of Fig. 33 is cut off at a value of 1.4 because the outlier frequencies for the highest-order approximations would make the remaining part of the plot completely unreadable. Fig. 35 shows the spectra obtained employing a uniformly spaced control net.

Vibrations of a clamped thin circular plate using three-dimensional solid elements
It was shown in Hughes et al. [27] that higher-order three-dimensional NURBS elements could be effectively utilized in the analysis of thin structures. In this section we consider the vibrations of a clamped, thin circular plate modeled as a three-dimensional solid. A coarse mesh, but one capable of exactly representing the geometry, is utilized and the order of the basis functions is increased by way of the k-refinement strategy (see [27]). The exact Poisson-Kirchhoff solution for this problem, given, for example, in [28],i s where R is the radius of the plate, t is the thickness, D ¼ Et 3 12ð1Àm 2 Þ is the flexural stiffness (E and m are Young's modulus and Poisson's ratio, resp.) and q is the density (mass per unit volume). For the first three frequencies, the values of the coefficients C mn are C 01 = 1.015,C 11 = 1.468 and C 02 = 2.007. The data for the problem are presented in Table 1. Note that, because the radius to thickness ratio is 100, the plate may be considered thin, and the results of Poisson-Kirchhoff theory may be considered valid.
The initial control net consists of 9 · 4 · 3 control points in the h, r,andz directions, respectively, and quadratic approximations in all the parametric directions are employed. Fig. 36 shows the mesh, consisting of eight elements within a single patch. The numerical results are compared with the exact solution in Table 2, where p, q and r are the orders of the basis functions in the circumferential, radial and vertical directions, respectively. Figs. 37-39 show the first three eigenmodes  (computed using p =4,q =5,r = 2), which are in qualitative agreement with the ones depicted in [28]. The relative errors (i.e., (x h À x)/x) for these cases are, respectively, 0.0054, 0.00027, and 0.0012.

NASA Aluminum Testbed Cylinder
Comparison studies between computational and experimental results are often humbling to analysts and experimentalists alike. Occasionally, large discrepancies occur. Sometimes, through careful reanalysis and/or careful re-experimentation, the root cause of the discrepancies can be identified but this is not always the case. Differences in ''as-built'' and ''as-drawn'' geometries, and ambiguities concerning material properties and boundary conditions, often make precise correlation impossible. In addition, errors can obviously be made in analysis. The same is true in experimentation. Nevertheless, comparison between calculations and experimental tests is fundamentally important in engineering. Here we compare frequencies calculated for exact, ''as-drawn'' models of structures for which considerable experimental data are available.   that control nets amount to a typical trilinear hexahedral mesh (see Figs. 48 and 53). This suggests the possibility that hexahedral finite element mesh generators (see, e.g., CUBIT [7]) may be useful in building isogeometric NURBS models. The stringer-main rib and stringer-end rib junctions are shown in Figs. 54 and 55, respectively. Note that there are gaps between the stringer and the ribs in the notch regions. Experimental vibration data have been obtained for a typical isolated stringer, a typical isolated main rib, the frame assembly, and the frame and skin assembly (see [23,9,18] for details and reference computational results). Calculations were performed of each of the corresponding isogeometric models. The Arnoldi Package (ARPACK, see [5]) eigensolver was used in the calculations of the individual components of the ATC,    Further refinement may be necessary to resolve the solution, as is the case with standard finite element analysis, but the geometry is never altered as the mesh is refined.  while the automated multi-level substructuring (AMLS) eigensolver (see [6]) was used in the calculations of the framework and the full ATC structure.
Frequency results for the stringer are presented in Fig. 56 and the first three bending modes are depicted in Fig. 57. The results shown are from analysis of a single patch with one rational quadratic element in the thickness, nine rational quadratic elements through the cross-section (the smallest number capable of exactly representing the geometry), and sixteen    C 5 rational sextic (p = 6) elements in the longitudinal direction. The resulting 144 element module has a total of 11,286 degrees-of-freedom. Other combinations of polynomial order, continuity and mesh size in the longitudinal direction were investigated. Higher-order, smooth meshes provided the most accuracy per degree-of-freedom but the overall run-time inevitably suffers if the polynomial order is raised indefinitely. Further study of the myriad of refinement options offered by NURBS-based isogeometric analysis is warranted.
The main rib frequency results are presented in Fig. 58. The ''coarse mesh'' (not shown) was comprised of 24 identical but rotated 15°sections, each comprised of six NURBS patches. Rational quadratic elements were used throughout. The full mesh for an individual rib had 34,704 degrees-of-freedom. Both h-andp-refinement were investigated for the reasons described below. In all cases, the results converged to the fine mesh results in Fig. 58.
The isolated main rib was the only case in which some of the numerically calculated frequencies were smaller than the experimental results. The fine mesh results fall below the coarse mesh results which, for theoretical reasons must occur (see,     e.g., [39]). The average error for the first eight modes is larger for the fine mesh than the coarse mesh, a somewhat surprising result. Due to the upper bound property of frequencies in our formulation, and the fact that our model is geometrically exact in the sense of the design drawings, we surmise there is some discrepancy between the drawings and the as-built configuration, or some other discrepancy between the experimental configuration and our model. Nevertheless, the correlation is still reasonable. Further study is needed to determine the cause of the differences. Selected modes shapes for the fine mesh are shown in Figs. 59 and 60.   Frequency results for the frame assembly are presented in Figs. 61 and 62. As for the case of the isolated stringer, the numerical results lie above the experimental results. The first bending and torsional modes of the frame assembly are shown in Figs. 63-66. A detail of the deformation pattern in the vicinity of a main rib-stringer junction for the first torsional mode   of the frame assembly is shown in Fig. 64. A mesh of 112,200 rational quadratic elements and 1,281,528 degrees-of-freedom was used for the analysis shown. One could reduce the number of degrees-of-freedom significantly by exploiting rotational symmetry and modeling only 1/24 of the frame assembly (as others have done, see [18]), but part of the goal of this work   was to demonstrate the feasibility of modeling an entire real structure of engineering interest using isoparametric NURBS elements, and so no such simplifications were employed.
Results for the frame and skin assembly are presented in Figs. 67 and 68. Once again, the numerical results lie above the experimental results. The first two modes are shown in Figs. 69-71. The mesh consisted of 228,936 rational quadratic elements and 2,219,184 degrees-of-freedom. The cost of array formation and assembly was commensurate with standard quadratic finite elements.

Conclusions
Isogeometric analysis was briefly reviewed and applied for the first time to structural vibrations. A number of elementary model problems were solved numerically and, in some cases, analytically. The models consisted of rods, beams, membranes, plates, and three-dimensional solids. Initiatory studies of rotationless bending elements were also presented.
The k-method was shown to provide more robust and accurate frequency spectra than typical higher-order finite elements (i.e., the p-method). Particularly intriguing is the possibility of eliminating ''optical'' branches of frequency spectra through the use of nonlinear parameterizations of the geometrical mapping. Optical branches have been identified as causes of severe accuracy degradation in higher modes and Gibbs phenomena in wave propagation. An exact geometrical model of the NASA Aluminum Testbed Cylinder was constructed and analyzed. Comparison was made with experimental results and good agreement was attained.
An important area of future research is mass lumping. A simple ''row sum'' procedure was investigated but it did not maintain the higher-order accuracy of consistent mass.

A.2. Order of accuracy employing cubic NURBS and consistent mass
Using cubic NURBS, the normalized discrete spectrum is given by Expanding cos(xh) and repeating the computations as before, we obtain: x h x $ 1 þ ðxhÞ 6 60; 480 ðA:9Þ and the order of accuracy in this case is 6. Fig. A.2 compares (A.9) with a spectrum obtained numerically.

A.3. Order of accuracy employing lumped mass
In similar fashion to the case of consistent mass, the order of accuracy of lumped mass can be obtained. Employing quadratic NURBS, we derive In these cases, the analytical expressions do not reproduce the behavior of the numerical spectra in the large, but only for the low-frequency part before the slope discontinuity, as shown in Fig. A.3. This is the relevant part as far as order of accuracy is concerned. So, by means of Taylor expansions, we obtain, for quadratic NURBS: x h x $ 1 À ðxhÞ 2 8 ðA:12Þ and for cubic NURBS x h x $ 1 À ðxhÞ 2 6 . ðA:13Þ As was already evident from  Rod problem and ''row sum'' lumped mass. Normalized discrete spectrum using cubic NURBS compared with 1 À (xh) 2 /6 for low frequencies.