# Isogeometric analysis of structural vibrations

J. Cottrell, R. Reali, Y. Bazilevs, T.J.R. Hughes

► To cite this version:

J. Cottrell, R. Reali, Y. Bazilevs, T.J.R. Hughes. Isogeometric analysis of structural vibrations. Computer Methods in Applied Mechanics and Engineering, 2006, 195, pp.5257-5296. 10.1016/j.cma.2005.09.027. hal-01516398

HAL Id: hal-01516398

https://hal.science/hal-01516398

Submitted on 30 Apr 2017

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Isogeometric analysis of structural vibrations

J.A. Cottrell, A. Reali, Y. Bazilevs, T.J.R. Hughes

Institute for Computational Engineering and Sciences, Computational and Applied Mathematics Chair III, The University of Texas at Austin,
201 East 24th Street, ACES 5.430A, 1 University Station C0200, Austin, TX 78712, USA

“Let no man ignorant of geometry enter,” – inscription above the entry to Plato’s Academy, circa 360 BC

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 fundamentals 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*-refinement is explored and shown to produce more accurate and robust results than corresponding finite elements. Through the use of nonlinear parameterization, “optical” branches of frequency spectra are eliminated for *k*-refined meshes. Optical branches have been identified as contributors to Gibbs phenomena in wave propagation problems and the cause of rapid degradation of higher modes in *p*-method finite 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.

**Keywords:** Isogeometric analysis; Structural vibrations; Invariant frequency spectra; Nonlinear parameterization; Finite elements; Rotationless bending elements; Mass matrices

## 1. 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
expressing the equation-by-equation duality of the force and displacement methods. John clearly had a flair for visual pre-
sentation of scientific ideas, and once stated that to do structural analysis well, one needed a sense of its “geometrical
beauty.” It is no wonder that John was initially attracted to structural analysis given that it is perhaps the most geometrical
of scientific subjects. This spirit is also vividly expressed in his book *Modern Fuselage Analysis and the Elastic Aircraft* [4].

I also remember a prescient article by Argyris entitled, *The Computer Shapes the Theory* [1]. It was clear John had precisely discerned the future whereas most others were oblivious, and remained so. It may come as a surprise to many young researchers today but, at that time, many scientific luminaries were fiercely resistant to the use of digital computers, and denigrated any work in which computers played an essential role. For example, the prominent mechanician, Clifford Truesdell, once gave a lecture in Milan entitled, “The Computer: Ruin of Science and Threat to Mankind” [40]. John had clearly foreseen and embraced the future, so attempts like this to legislate science must have seemed particularly laughable.

John's early work attracted me to the subject but I never met him until the first FENOMECH conference at the Uni-
versity of Stuttgart in 1978 [36]. He was very gregarious and friendly when I met him. It was not very long after our first
meeting that he invited me to join him and William Prager as co-editor of this journal. Frankly, I was shocked at the time
because it was only a few years after I completed my Ph.D. Before beginning my duties as co-editor, William Prager passed
away. Thereafter, John, Tinsley Oden and I worked together as the editorial team until John retired. During those years,
my relationship with John was always very friendly. We spoke on the telephone about editorial issues and had meetings
together with the publisher at scientific conferences. He was always respectful of my ideas and gave me complete autonomy
over my editorial activities. I always received Christmas cards from him with personal notes and an occasional letter. A
frequent theme was the celebration of our friendship. I have fond memories of him.

As I think back about the early days of finite elements, and John Argyris's work in particular, I am struck by the empha-
sis on geometry. This is understandable because geometry is the foundation of analysis and, at the time, there was no com-
puter representation of geometry that could support computer based analysis. Geometry was encapsulated in mechanical
drawings and these were done with pencils and pens on paper, Mylar and vellum.<fnref n="1" /> Consequently, computer representation
of geometry had to be invented and the finite element description was, first and foremost, simply a geometry. With inspi-
ration from the matrix theory of frame structures, the idea of geometric discretization of two- and three-dimensional struc-
tures into “elements” was in the air. In fact, this concept seems to have taken hold before there was general agreement as to
how to derive element stiffness matrices for even the simplest elements. It is interesting to recall the premise behind the 1964
monograph [22]. Its author, Richard H. Gallagher, another finite element pioneer, set out to determine how stiffness matri-
ces for the four-node parallelogram, plane stress element were being calculated at various centers around the world
engaged in finite element research. Apparently, no two results were identical. Everyone understood that basic element
shapes were necessary to build the computer-based geometrical descriptions. The analytical underpinnings followed.
(Indeed, “The computer shapes the theory!”)

Computer aided design (CAD) came later, in the 1970s and 1980s, and with it came a new and different geometrical description of a structure. CAD and the finite element method have had a difficult relationship. They have never been married, and attempts at living together have also not been very successful. There seem to be irreconcilable differences in their geometric descriptions. CAD surface geometries frequently have gaps and overlaps which need to be fixed prior to finite element mesh generation. Often, features, such as holes, rivets, and bolts, need to be removed in order to generate meshes. Sometimes, features, such as weld fillets, have to be added in order to avoid physically unrealistic stresses. One might conclude that the CAD description of geometry is not a convenient basis for finite element analysis but one might also conclude that the finite element description of geometry, typically consisting of faceted, low-order polynomial elements, is not an adequate description of geometry for modern analysis. It is hard to imagine that analysis, on its own, could re-shape CAD because the latter is a significantly larger business than the former. However, CAD descriptions are also unsuitable for animation and visualization, and the computer gaming industry is significantly larger than CAD. So perhaps some changes may be possible in the ensuing years. Engineering analysis has not yet been impacted by major advances in computational geometry but the authors believe this is about to change because the time has come for a better geometrical description than that provided by finite elements. This paper takes a step in this direction. It is inspired by, and dedicated to, the memory of John Argyris.

**2. 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 func-
tions employed in computer aided design (CAD) systems, capable of representing many engineering geometries exactly; to

<footnote n="1">¹ Young engineers may not know what vellum and Mylar are. Vellum is a translucent drafting material made from cotton fiber. It is strong, erasable, and provides archival quality. Mylar is the trade name of a translucent polyester film used for drafting.</footnote>
<figure><img src="image_2.png" /><figcaption>Fig. 1. Schematic illustration comparing finite element analysis and isogeometric analysis meshes for a bracket.</figcaption></figure>

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 mod-
eled as three-dimensional solids, and particularly intriguing results were obtained for the linear advection–diffusion equa-
tion. 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–16,19,31–34]).<fnref n="2" />

<footnote n="2">² It is interesting to note that the BOSOR (Buckling of Shells of Revolution) codes, which were developed in the 1960s by David Bushnell at the Lockheed Palo Alto Research Laboratory, and still enjoy wide use, also employed a rotationless formulation (see [10,11]).</footnote>
Finite Element Analysis

<figure><img src="image_3.png" /></figure>

Isogeometric Analysis

<figure><img src="image_4.png" /><figcaption>Fig. 2. “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.</figcaption></figure>

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* three-dimensional 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.

### 3. 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].
### 3.1. 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.

#### 3.1.1. Knot vectors

A knot vector, $\Xi$, is a set of non-decreasing real numbers representing coordinates in the parametric space of the curve:

$$ \Xi = \{\xi_1, \dots, \xi_{n+p+1}\}, \qquad (1) $$

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 $[\xi_1, \xi_{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 $[\xi_1, \xi_{n+p+1}]$ but are not, in general, interpolatory at interior knots.

#### 3.1.2. Basis functions

Given a knot vector, $\Xi$, B-Spline basis functions are defined recursively starting with $p=0$ (piecewise constants):

$$ N_{i,0}(\xi) = \begin{cases} 1 & \text{if } \xi_i \le \xi < \xi_{i+1}, \\ 0 & \text{otherwise.} \end{cases} \qquad (2) $$

For $p > 1$:

$$ N_{i,p}(\xi) = \frac{\xi - \xi_i}{\xi_{i+p} - \xi_i} N_{i,p-1}(\xi) + \frac{\xi_{i+p+1} - \xi}{\xi_{i+p+1} - \xi_{i+1}} N_{i+1,p-1}(\xi). \qquad (3) $$

In Fig. 3 we present an example consisting of $n=9$ cubic basis functions generated from the open knot vector $\Xi = \{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, $\sum_{i=1}^{n} N_{i,p}(\xi) = 1 \quad \forall \xi$.

* The support of each $N_{i,p}$ is compact and contained in the interval $[\xi_i, \xi_{i+p+1}]$.

* B-Spline basis functions are non-negative: $N_{i,p}(\xi) \ge 0 \quad \forall \xi$.

#### 3.1.3. 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 $\mathbb{R}^d$, the piecewise polynomial B-Spline curve $C(\xi)$ of order $p$ is obtained by taking a linear combination of basis functions and control points:

$$ C(\xi) = \sum_{i=1}^{n} N_{i,p}(\xi) B_i, \qquad (4) $$

<figure><img src="image_5.png" /><figcaption>Fig. 3. Cubic basis functions formed from the open knot vector $\Xi = \{0, 0, 0, 0, 1/6, 1/3, 1/2, 2/3, 5/6, 1, 1, 1, 1\}$.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>Fig. 4. Piecewise cubic B-Spline curve (solid line) and its control net (dotted).</figcaption></figure>

where $\mathbf{B}_i$ is the $i$th control point. The piecewise linear interpolation of the control points defines the control net. In Fig. 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.

### 3.1.4. B-Spline surfaces

By means of tensor products, B-Spline surfaces can be constructed starting from knot vectors $\Xi = \{\xi_1, \dots, \xi_{n+p+1}\}$ and $H = \{\eta_1, \dots, \eta_{m+q+1}\}$, and an $n \times m$ net of control points $\mathbf{B}_{i,j}$. One-dimensional basis functions $N_{i,p}$ and $M_{j,q}$ (with $i=1, \dots, n$ and $j=1, \dots, m$) of order $p$ and $q$, respectively, are defined from the knot vectors, and the B-Spline surface is constructed as

$$S(\xi, \eta) = \sum_{i=1}^{n} \sum_{j=1}^{m} N_{i,p}(\xi) M_{j,q}(\eta) \mathbf{B}_{i,j}. \qquad (5)$$

The patch is the domain $[\xi_1, \xi_{n+p+1}] \times [\eta_1, \eta_{m+q+1}]$.

### 3.1.5. B-Spline solids

B-Spline solids can be constructed in similar fashion. Given an $n \times m \times l$ control net and three knot vectors:

$$\Xi = \{\xi_1, \dots, \xi_{n+p+1}\}, \quad H = \{\eta_1, \dots, \eta_{m+q+1}\}, \quad Z = \{\zeta_1, \dots, \zeta_{l+r+1}\} \qquad (6)$$

from which the one-dimensional basis functions $N_{i,p}$, $M_{j,q}$ and $L_{k,r}$ (with $i=1, \dots, n$, $j=1, \dots, m$ and $k=1, \dots, l$) of order $p$, $q$ and $r$, respectively, are defined, the B-Spline solid is then:

$$S(\xi, \eta, \zeta) = \sum_{i=1}^{n} \sum_{j=1}^{m} \sum_{k=1}^{l} N_{i,p}(\xi) M_{j,q}(\eta) L_{k,l}(\zeta) \mathbf{B}_{i,j,k}. \qquad (7)$$

## 3.2. Non-uniform rational B-Splines

A rational B-Spline in $\mathbb{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 $\mathbb{R}^d$, we start from a set $\mathbf{B}_i^w$ ($i=1, \dots, n$) of control points (“projective points”) for a B-Spline curve in $\mathbb{R}^{d+1}$ with knot vector $\Xi$. Then the control points for the NURBS curve are

$$(\mathbf{B}_i)_j = \frac{(\mathbf{B}_i^w)_j}{w_i}, \quad j = 1, \dots, d, \qquad (8)$$
where $(\mathbf{B}_i)_j$ is the jth component of the vector $\mathbf{B}_i$ and $w_i = (\mathbf{B}_i^w)_{d+1}$ is referred to as the ith weight. The NURBS basis functions of order p are then defined as

$$R_i^p(\xi) = \frac{N_{i,p}(\xi)w_i}{\sum_{i=1}^n N_{i,p}(\xi)w_i} \qquad (9)$$

and their first and second derivatives are given by

$$(R_i^p)'(\xi) = \frac{N_{i,p}'(\xi)w_i}{\sum_{i=1}^n N_{i,p}(\xi)w_i} - \frac{N_{i,p}(\xi)w_i \sum_{i=1}^n N_{i,p}'(\xi)w_i}{\left(\sum_{i=1}^n N_{i,p}(\xi)w_i\right)^2} \qquad (10)$$

and

$$(R_i^p)''(\xi) = \frac{N_{i,p}''(\xi)w_i}{\sum_{i=1}^n N_{i,p}(\xi)w_i} + \frac{2N_{i,p}(\xi)w_i \left(\sum_{i=1}^n N_{i,p}'(\xi)w_i\right)^2}{\left(\sum_{i=1}^n N_{i,p}(\xi)w_i\right)^3} - \frac{2N_{i,p}'(\xi)w_i \sum_{i=1}^n N_{i,p}(\xi)w_i + N_{i,p}(\xi)w_i \sum_{i=1}^n N_{i,p}''(\xi)w_i}{\left(\sum_{i=1}^n N_{i,p}(\xi)w_i\right)^2} \qquad (11)$$

The NURBS curve is defined by

$$\mathbf{C}(\xi) = \sum_{i=1}^{n} R_{i}^{p}(\xi)\mathbf{B}_{i}. \qquad (12)$$

Rational surfaces and solids are defined in an analogous way in terms of the basis functions (resp.)

$$R_{i,j}^{p,q}(\xi, \eta) = \frac{N_{i,p}(\xi)M_{j,q}(\eta)w_{i,j}}{\sum_{i=1}^{n} \sum_{j=1}^{m} N_{i,p}(\xi)M_{j,q}(\eta)w_{i,j}} \qquad (13)$$

and

$$R_{i,j,k}^{p,q,r}(\xi, \eta, \zeta) = \frac{N_{i,p}(\xi)M_{j,q}(\eta)L_{k,r}(\zeta)w_{i,j,k}}{\sum_{i=1}^{n} \sum_{j=1}^{m} \sum_{k=1}^{l} N_{i,p}(\xi)M_{j,q}(\eta)L_{k,r}(\zeta)w_{i,j,k}} \qquad (14)$$

In the following, we summarize noteworthy properties of NURBS:

* NURBS basis functions formed from an open knot vector constitute a partition of unity: $\sum_{i=1}^n R_i^p(\xi) = 1 \forall \xi$.

* 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.

### 3.3. 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 $\Xi \times H \times 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 p-refinement (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).

### 3.4. *k*-Refinement

Isogeometric analysis is fundamentally a higher-order approach. While it is true that the first two NURBS bases consist of constants and linear, 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*- and *p*-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 *k*-refinement, 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-

<figure><img src="image_7.png" /><figcaption>Fig. 5. Comparison of control variable growth in one dimension.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Fig. 6. Comparison of control variable growth in two dimensions.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>Fig. 7. Comparison of control variable growth in three dimensions.</figcaption></figure>

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 d-dimensional 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 + dr^{n-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.

# 4. Structural vibrations

We briefly recall the main equations of structural vibrations; for elaboration, see [26,17,12].

## 4.1. Natural frequencies and modes

Given a linear multi-degree-of-freedom structural system, the undamped, unforced equations of motion, which govern free vibrations, are

$$M\ddot{u} + Ku = 0, \qquad (15)$$

where M and K are, respectively, the consistent mass and the stiffness matrices, $u = u(t)$ is the displacement vector and $\ddot{u} = \frac{d^2u}{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, $\phi_n$, can be obtained by separation of variables:

$$u(t) = \phi_n e^{i\omega_n t}, \qquad (16)$$

where $\omega_n$ is the nth natural frequency. Combining (15) and (16) leads to the generalized eigenproblem:

$$(K - \omega_n^2 M)\phi_n = 0. \qquad (17)$$

The normal modes are defined up to a multiplicative constant. Different ways of normalization have been proposed. One of the most widely used is

$$\phi_n^T M \phi_n = 1. \qquad (18)$$

# 5. 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).
## 5.1. 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

$$ 
\begin{aligned}
u_{xx} + \omega^2 u &= 0 \quad \text{for } x \in ]0, 1[ \\
u(0) &= u(1) = 0
\end{aligned} 
\qquad (19) 
$$

and for which the exact solution in terms of natural frequencies is

$$ \omega_n = n\pi \quad \text{with } n = 1, 2, 3, \dots \qquad (20) $$

After writing the weak formulation and performing the discretization, a problem of the form of (17) is obtained.

### 5.1.1. 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, $\omega_n^h$, 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 degrees-of-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$. In Fig. 10, this behavior is highlighted by plotting the last computed frequencies for $p = 2, \dots, 10$. To understand the outliers, we first remark that the finite element spectrum for quadratic elements consists of acoustic (low-mode) and optical (high-mode) 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–26]). However, they would be detrimental in explicit transient analysis because the frequencies of

<figure><img src="image_10.png" /><figcaption>Fig. 8. Fixed-fixed rod. Normalized discrete spectra using quadratic finite elements and NURBS.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>Fig. 9. Fixed-fixed rod. Normalized discrete spectra using different order NURBS basis functions. Outliers appear in the very thin band on the right end of the spectrum.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Fig. 10. Fixed-fixed rod. Last normalized frequencies for $p = 2, \dots, 10$.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Fig. 11. Fixed-fixed rod. Order of convergence for the first three frequencies using quadratic NURBS.</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>Fig. 12. Fixed-fixed rod. Order of convergence for the first three frequencies using cubic NURBS.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Fig. 13. Fixed-fixed rod. Order of convergence for the first three frequencies using quartic NURBS.</figcaption></figure>

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.

### 5.1.2. 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

$$ \mathbf{M}^e = \frac{h}{120} \begin{bmatrix} 6 & 13 & 1 \\ 13 & 54 & 13 \\ 1 & 13 & 6 \end{bmatrix}, \quad \mathbf{K}^e = \frac{1}{6h} \begin{bmatrix} 2 & -1 & -1 \\ -1 & 2 & -1 \\ -1 & -1 & 2 \end{bmatrix}, \qquad (21) $$

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$, is

$$ \frac{h}{120} (\ddot{u}_{A-2} + 26\ddot{u}_{A-1} + 66\ddot{u}_{A} + 26\ddot{u}_{A+1} + \ddot{u}_{A+2}) - \frac{1}{6h} (u_{A-2} + 2u_{A-1} - 6u_{A} + 2u_{A+1} + u_{A+2}) = 0, \qquad (22) $$
which can be compactly written as

$$ \frac{h^2}{20} \alpha \ddot{u}_A - \beta u_A = 0, \qquad (23) $$

where $\alpha$ and $\beta$ are operators defined as follows:

$$ \begin{aligned} \alpha x_A &= x_{A-2} + 26x_{A-1} + 66x_A + 26x_{A+1} + x_{A+2}, \\ \beta x_A &= x_{A-2} + 2x_{A-1} - 6x_A + 2x_{A+1} + x_{A+2}. \end{aligned} \qquad (24) $$

Separating variables

$$ u_A(t) = \phi_A q(t) \qquad (25) $$

and substituting this expression into (23), after adding and subtracting $\frac{(\omega^h h)^2}{20} \alpha u_A$, we obtain:

$$ \left( \ddot{q} + (\omega^h)^2 q \right) \frac{h^2}{20} \alpha \phi_A - \frac{(\omega^h h)^2}{20} \alpha \phi_A + \beta \phi_A \right) q = 0. \qquad (26) $$

The satisfaction of (26) is achieved by selecting $\phi_A$ and $q$ such that:

$$ \frac{(\omega^h h)^2}{20} \alpha + \beta \phi_A = 0 \qquad (27) $$

and

$$ \ddot{q} + (\omega^h)^2 q = 0. \qquad (28) $$

Assuming a solution for (27) of the form (for fixed-fixed boundary conditions):

$$ \phi_A = C \sin(A\omega h), \quad \omega = n\pi, \qquad (29) $$

Eq. (27) can be rewritten as

$$ \frac{(\omega^h h)^2}{20} \alpha + \beta \sin(A\omega h) = 0. \qquad (30) $$

Substituting expressions (24) for $\alpha$ and $\beta$, and using the trigonometric identity $\sin(a \pm b) = \sin(a) \cos(b) \pm \sin(b) \cos(a)$, yields:

$$ \frac{(\omega^h h)^2}{20} (16 + 13 \cos(\omega h) + \cos^2(\omega h)) - (2 - \cos(\omega h) - \cos^2(\omega h)) = 0, \qquad (31) $$

which can be solved for $\frac{\omega^h}{\omega}$, giving

$$ \frac{\omega^h}{\omega} = \frac{1}{\omega h} \sqrt{\frac{20(2 - \cos(\omega h) - \cos^2(\omega h))}{16 + 13 \cos(\omega h) + \cos^2(\omega h)}}. \qquad (32) $$

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

$$ \frac{\omega^h}{\omega} = \frac{1}{\omega h} \sqrt{\frac{42(16 - 3 \cos(\omega h) - 12 \cos^2(\omega h) - \cos^3(\omega h))}{272 + 297 \cos(\omega h) + 60 \cos^2(\omega h) + \cos^3(\omega h)}}. \qquad (33) $$

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 1.** Eqs. (32) and (33), and Fig. 14, confirm that the continuous part of the NURBS frequency spectra are invariant, that is, are independent of the number of degrees-of-freedom.

**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
<figure><img src="image_16.png" /><figcaption>Fig. 14. Fixed-fixed rod. Comparison of analytical and numerical discrete spectra computed for quadratic and cubic NURBS.</figcaption></figure>

presented in Fig. 15. This choice corresponds to a nonlinear parameterization (see Figs. 16 and 17 for plots of the parameterization $x(\xi)$ and its Jacobian $J(\xi) = \frac{dx(\xi)}{d\xi}$ 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 for quadratic and cubic NURBS.

Outliers, optical branches, and nonlinear parameterizations are curious entities. Their practical implications are not well understood at this time.

<figure><img src="image_17.png" /><figcaption>Fig. 15. Control points for linear parameterization (dots) compared with uniformly spaced control points (asterisks) for cubic NURBS.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>Fig. 16. Plot of the parameterization for the cases of uniformly spaced control points and linear parameterization (cubic NURBS, 21 control points).</figcaption></figure>
<figure><img src="image_19.png" /><figcaption>Fig. 17. Plot of the Jacobian of the parameterization for the cases of uniformly spaced control points and linear parameterization (cubic NURBS, 21 control points).</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>Fig. 18. Fixed-fixed rod. Normalized discrete spectra using uniformly spaced control points. (These results are identical to those presented in Fig. 9 except the outliers have been eliminated.)</figcaption></figure>

**Remark 3.** In this work, consistent mass is emphasized because it seems more suitable than lumped mass when higher-order 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.

## 5.2. 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

$$ 
\begin{aligned}
u_{xxxx} - \omega^2 u = 0 \quad & \text{for } x \in ]0, 1[, \\
u(0) = u(1) = u_{xx}(0) = u_{xx}(1) = 0,
\end{aligned} 
\qquad (34) $$
<figure><img src="image_21.png" /><figcaption>Fig. 19. Fixed-fixed rod. Normalized discrete spectra using different order NURBS basis functions with “row sum” lumped mass matrices.</figcaption></figure>

where

$$ \omega_n = (n\pi)^2 \quad \text{with } n = 1, 2, 3, \dots \qquad (35) $$

### 5.2.1. 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–25 show that the order of convergence of frequencies using NURBS is optimal, that is $O(h^{2(p-1)})$.

<figure><img src="image_22.png" /><figcaption>Fig. 20. Simply supported beam. Normalized discrete spectra using cubic finite elements and NURBS.</figcaption></figure>
<figure><img src="image_23.png" /><figcaption>Fig. 21. Simply supported beam. Normalized discrete spectra using different order NURBS basis functions.</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>Fig. 22. Simply supported beam. Last normalized frequencies for $p = 2, \dots, 10$.</figcaption></figure>

<figure><img src="image_25.png" /><figcaption>Fig. 23. Simply supported beam. Order of convergence for the first three frequencies using quadratic NURBS.</figcaption></figure>
<figure><img src="image_26.png" /><figcaption>Fig. 24. Simply supported beam. Order of convergence for the first three frequencies using cubic NURBS.</figcaption></figure>

<figure><img src="image_27.png" /><figcaption>Fig. 25. Simply supported beam. Order of convergence for the first three frequencies using quartic NURBS.</figcaption></figure>

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:

$$ \frac{\omega^h}{\omega} = \frac{1}{\omega h^2} \sqrt{\frac{210(2 - 3 \cos(\omega h) + \cos^3(\omega h))}{272 + 297 \cos(\omega h) + 60 \cos^2(\omega h) + \cos^3(\omega h)}} \quad (36) $$

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.

### 5.2.2. 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:
<figure><img src="image_28.png" /><figcaption>Fig. 26. Simply supported beam. Analytical and numerical discrete spectra computed using cubic and quartic NURBS.</figcaption></figure>

<figure><img src="image_29.png" /><figcaption>Fig. 27. Simply supported beam. Normalized discrete spectra using equally spaced control points. (These results are identical to those presented in Fig. 21 except the outliers are eliminated.)</figcaption></figure>

$$ \begin{aligned} u_{xxxx} - \omega^2 u = 0 \quad & \text{for } x \in [0, 1[, \\ u(0) = u_x(0) = u_{xx}(1) = u_{xxx}(1) = 0. \end{aligned} \tag{37} $$

The natural frequencies are (see [12]) $\omega_n = \beta_n^2$, with $\beta_1 = 1.8751$, $\beta_2 = 4.6941$, $\beta_3 = 7.8548$, $\beta_4 = 10.996$, and $\beta_n = (n - 1/2)\pi$ 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

$$ A(v^h, u^h) = \int_0^1 v_{xx}^h u_{xx}^h \, dx + v_x^h u_{xx}^h|_{x=0} + v_{xx}^h u_x^h|_{x=0} + \tau v_x^h u_x^h|_{x=0}, \tag{38} $$

where $v^h$ and $u^h$ are the discrete weighting and trial solution, respectively, and $\tau$ 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 $\tau$ 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 $\tau = 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

$$ A(v^h, \mu, u^h, \lambda) = \int_0^1 v_{xx}^h u_{xx}^h \, dx + \lambda v_x^h|_{x=0} + \mu u_x^h|_{x=0}, \tag{39} $$
<figure><img src="image_30.png" /><figcaption>Fig. 28. Cantilever beam with weakly enforced rotation boundary condition. Normalized discrete spectra using different order NURBS basis functions.</figcaption></figure>

<figure><img src="image_31.png" /><figcaption>Fig. 29. Cantilever beam with Lagrange multiplier. Normalized discrete spectra using different order NURBS basis functions.</figcaption></figure>

where $\lambda$ is the Lagrange multiplier and $\mu$ 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).

# 6. 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.

## 6.1. 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

$$ 
\begin{aligned}
\nabla^2 u(x, y) + \omega^2 u(x, y) &= 0, \quad (x, y) \in \Omega =]0, 1[\times]0, 1[, \\
u(x, y)|_{\partial\Omega} &= 0,
\end{aligned} 
\quad (40) $$
where $\nabla^2$ is the Laplace operator. The exact natural frequencies are (see, e.g., [28]):

$$ \omega_{mn} = \pi\sqrt{m^2 + n^2}, \quad m, n = 1, 2, 3, \dots \qquad (41) $$

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.

<figure><img src="image_32.png" /><figcaption>Fig. 30. Square membrane. Normalized discrete spectra using different order NURBS basis functions (90 × 90 control points). Note the presence of optical branches.</figcaption></figure>

<figure><img src="image_33.png" /><figcaption>Fig. 31. Square membrane. Detail of the low-frequency part of the normalized discrete spectra.</figcaption></figure>
<figure><img src="image_34.png" /><figcaption>Fig. 32. Square membrane. Normalized discrete spectra using a uniformly spaced control net. Note, the optical branches of Fig. 30 are eliminated.</figcaption></figure>

## 6.2. 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:

$$ \nabla^4 u(x,y) - \omega^2 u(x,y) = 0 \quad \text{for } (x,y) \in \Omega = ]0, 1[ \times ]0, 1[, \quad (42) $$

$$ u(x,y)|_{\partial\Omega} = 0, \quad \nabla^2 u(x,y)|_{\partial\Omega} = 0, $$

for which the exact solution natural frequencies (see, e.g., [28]) are

$$ \omega_{mn} = \pi^2(m^2 + n^2), \quad m, n = 1, 2, 3, \dots \quad (43) $$

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).

<figure><img src="image_35.png" /><figcaption>Fig. 33. Poisson–Kirchhoff plate. Normalized discrete spectra using different order NURBS basis functions (90 × 90 control points). Note the presence of optical branches.</figcaption></figure>
<figure><img src="image_36.png" /><figcaption>Fig. 34. Poisson–Kirchhoff plate. Detail of the low-frequency part of the normalized discrete spectra.</figcaption></figure>

<figure><img src="image_37.png" /><figcaption>Fig. 35. Poisson–Kirchhoff plate. Normalized discrete spectra using a uniformly spaced control net. Note, the optical branches of Fig. 33 are eliminated.</figcaption></figure>

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.

## 7. 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], is
<table><caption>Table 1
Clamped circular plate, geometric and material parameters</caption><tr><td>R</td><td>2 m</td></tr><tr><td>t</td><td>0.02 m</td></tr><tr><td>E</td><td>30 × 10<sup>6</sup> KN/m<sup>2</sup></td></tr><tr><td>ν</td><td>0.2</td></tr><tr><td>ρ</td><td>2.320 KN s/m<sup>4</sup></td></tr></table>

$$ \omega_{mn} = C_{mn}^2 \frac{\pi^2}{R^2} \sqrt{\frac{D}{\rho t}} \quad [\text{rad/s}], \qquad (44) $$

where R is the radius of the plate, t is the thickness, $D = \frac{Et^3}{12(1-ν^2)}$ is the flexural stiffness (E and ν are Young's modulus and Poisson's ratio, resp.) and ρ 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 θ, r, and z 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

<figure><img src="image_38.png" /><figcaption>Fig. 36. Clamped circular plate. Eight element mesh.</figcaption></figure>

<table><caption>Table 2
Clamped circular plate, numerical results compared with the exact solution</caption><thead><tr><th>p</th><th>q</th><th>r</th><th>ω<sub>01</sub> [rad/s]</th><th>ω<sub>11</sub> [rad/s]</th><th>ω<sub>02</sub> [rad/s]</th></tr></thead><tbody><tr><td>2</td><td>2</td><td>2</td><td>138.133</td><td>1648.800</td><td>2052.440</td></tr><tr><td>2</td><td>3</td><td>2</td><td>56.702</td><td>267.765</td><td>276.684</td></tr><tr><td>3</td><td>3</td><td>2</td><td>56.051</td><td>126.684</td><td>232.788</td></tr><tr><td>3</td><td>4</td><td>2</td><td>54.284</td><td>124.417</td><td>212.451</td></tr><tr><td>4</td><td>4</td><td>2</td><td>54.284</td><td>113.209</td><td>212.451</td></tr><tr><td>4</td><td>5</td><td>2</td><td>54.153</td><td>112.700</td><td>210.840</td></tr><tr><td>Exact</td><td></td><td></td><td>53.863</td><td>112.670</td><td>210.597</td></tr></tbody></table>

<figure><img src="image_39.png" /><figcaption>Fig. 37. Clamped circular plate. Eigenmode corresponding to $ω_{01}$.</figcaption></figure>
<figure><img src="image_40.png" /><figcaption>Fig. 38. Clamped circular plate. Eigenmode corresponding to $\omega_{11}$.</figcaption></figure>

<figure><img src="image_41.png" /><figcaption>Fig. 39. Clamped circular plate. Eigenmode corresponding to $\omega_{02}$.</figcaption></figure>

(computed using $p = 4, q = 5, r = 2$), which are in qualitative agreement with the ones depicted in [28]. The relative errors (i.e., $(\omega^h - \omega)/\omega$) for these cases are, respectively, 0.0054, 0.00027, and 0.0012.

## 8. 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.

<figure><img src="image_42.png" /><figcaption>Fig. 40. NASA Aluminum Testbed Cylinder (ATC). Frame and skin.</figcaption></figure>
The NASA Aluminum Testbed Cylinder (ATC) is shown in Figs. 40 and 41. An isogeometric model (see Fig. 42) was constructed from design drawings. There are three distinct members composing the framework (see Fig. 43): nine identical main ribs (see Figs. 44–48); 24 identical, prismatic stringers (see Fig. 49); and two end ribs, which are mirror images of each other (see Figs. 50–53). Note that every geometrical feature is *exactly* represented in the model. It is also interesting to note

<figure><img src="image_43.png" /><figcaption>Fig. 41. NASA ATC. Frame only.</figcaption></figure>

<figure><img src="image_44.png" /><figcaption>Fig. 42. NASA ATC frame and skin: isogeometric model.</figcaption></figure>

<figure><img src="image_45.png" /><figcaption>Fig. 43. NASA ATC frame: isogeometric model.</figcaption></figure>
<figure><img src="image_46.png" /><figcaption>Fig. 44. NASA ATC. Isogeometric model of the main rib.</figcaption></figure>

<figure><img src="image_47.png" /><figcaption>Fig. 45. NASA ATC. Typical 15° segment of the main rib. Mesh 1, the coarsest mesh, encapsulates the exact geometry.</figcaption></figure>

<figure><img src="image_48.png" /><figcaption>Fig. 46. NASA ATC. Typical 15° segment of the main rib. Mesh 2. Knot insertion has been used selectively to help even out the aspect ratios of elements making Mesh 2 more uniform and suitable for analysis.</figcaption></figure>

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,
<figure><img src="image_49.png" /><figcaption>Fig. 47. NASA ATC. Typical 15° segment of the main rib. Mesh 3. 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.</figcaption></figure>

<figure><img src="image_50.png" /><figcaption>Fig. 48. NASA ATC. Detail of the “notch” region in the main rib. The control net is on the left and the exact geometry is on the right.</figcaption></figure>

<figure><img src="image_51.png" /><figcaption>Fig. 49. NASA ATC. Isogeometric model of the longitudinal stringer. Sample meshes.</figcaption></figure>

<figure><img src="image_52.png" /><figcaption>Fig. 50. NASA ATC. Typical 15° segment of an end rib. Mesh 1 (coarsest mesh).</figcaption></figure>
<figure><img src="image_53.png" /><figcaption>Fig. 51. NASA ATC. Typical 15° segment of an end rib. Mesh 2.</figcaption></figure>

<figure><img src="image_54.png" /><figcaption>Fig. 52. NASA ATC. Typical 15° segment of an end rib. Mesh 3.</figcaption></figure>

<figure><img src="image_55.png" /><figcaption>Fig. 53. NASA ATC. Detail of the “notch” region in an end rib. The control net is on the left and the exact geometry is on the right.</figcaption></figure>

<figure><img src="image_56.png" /><figcaption>Fig. 54. NASA ATC. Stringer–main rib junction.</figcaption></figure>

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
<figure><img src="image_57.png" /><figcaption>Fig. 55. NASA ATC. Stringer-end rib junction.</figcaption></figure>

<figure><img src="image_58.png" /><figcaption>Fig. 56. NASA ATC. Comparison of numerical and experimental frequency results for the longitudinal stringer.</figcaption></figure>

<figure><img src="image_59.png" /><figcaption>Fig. 57. NASA ATC. Selected calculated mode shapes for the stringer. Three lowest x-z modes.</figcaption></figure>

$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*- and *p*-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,
<figure><img src="image_60.png" /><figcaption>Fig. 58. NASA ATC. Comparison of numerical and experimental frequency results for the main rib.</figcaption></figure>

<figure><img src="image_61.png" /><figcaption>Fig. 59. NASA ATC. Computed mode shapes for the main rib. First three out-of-plane modes.</figcaption></figure>
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

<figure><img src="image_62.png" /><figcaption>Fig. 60. NASA ATC. Computed mode shapes for the main rib. First three in-plane modes.</figcaption></figure>

<figure><img src="image_63.png" /><figcaption>Fig. 61. NASA ATC. Comparison of numerical and experimental frequency results for the frame assembly.</figcaption></figure>
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.

<figure><img src="image_64.png" /><figcaption>Fig. 62. NASA ATC. Relative frequency error for the frame assembly.</figcaption></figure>

<figure><img src="image_65.png" /><figcaption>Fig. 63. NASA ATC. Calculated first torsional mode for the frame assembly; side view. The color contours represent the vertical displacement.</figcaption></figure>

<figure><img src="image_66.png" /><figcaption>Fig. 64. NASA ATC. Detail of first torsional mode for the frame assembly; stringer–main rib junction.</figcaption></figure>
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

<figure><img src="image_67.png" /><figcaption>Fig. 65. NASA ATC. Calculated first torsional mode for the frame assembly; end view. The color contours represent the vertical displacement.</figcaption></figure>

<figure><img src="image_68.png" /><figcaption>Fig. 66. NASA ATC. Calculated first bending mode for the frame assembly. The color contours represent the vertical displacement.</figcaption></figure>

<figure><img src="image_69.png" /><figcaption>Fig. 67. NASA ATC. Comparison of numerical and experimental frequency results for the frame and skin assembly.</figcaption></figure>
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

<figure><img src="image_70.png" /><figcaption>Fig. 68. NASA ATC. Relative frequency error for the frame skin assembly.</figcaption></figure>

<figure><img src="image_71.png" /><figcaption>Fig. 69. NASA ATC. Calculated first Rayleigh mode of the frame and skin assembly. The color contours represent the ovalization of the assembly.</figcaption></figure>

<figure><img src="image_72.png" /><figcaption>Fig. 70. NASA ATC. Calculated first Love mode of the frame and skin assembly. The color contours represent the ovalization of the assembly.</figcaption></figure>
<figure><img src="image_73.png" /><figcaption>Fig. 71. NASA ATC. Calculated first Love mode of the frame and skin assembly. The color contours represent the axial displacement of the assembly.</figcaption></figure>

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.

## 9. 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.

## Acknowledgments

The authors wish to thank Jeff Bennighof and his research assistants for providing us with access to the AMLS eigensolver, which was used to calculate the natural frequencies and modes of the ATC. In particular, Mark Muller gave generously of his time when we were facing deadlines. We also wish to thank Ted Belytchko for pointing out the Brillouin reference to us, and Saikat Dey for providing us with the ATC data.

Alessandro Reali was a visiting scholar at the Institute for Computational Engineering and Sciences where he performed some of the research reported herein in partial fulfillment of the Masters Degree in Earthquake Engineering at the ROSE School (European School for Advanced Studies in Reduction of Seismic Risk), Istituto Universitario di Studi Superiori di Pavia, Università degli Studi di Pavia, Italy. J. Austin Cottrell was supported by a VIGRE Fellowship. Funding for this research was provided by the Office of Naval Research Contract N00014-03-0263, Dr. Luise Couchman, Contract Monitor. This support is gratefully acknowledged.

## Appendix A. Computation of the order of accuracy for the rod problem

Starting from the analytical expressions for the normalized discrete spectra obtained previously, it is possible to compute the corresponding orders of accuracy by means of Taylor expansions. The computation of the orders of accuracy for the rod problem using both consistent and lumped mass formulations and employing quadratic and cubic NURBS is presented as follows.
A.1. Order of accuracy employing quadratic NURBS and consistent mass

The analytical expression for the normalized discrete spectrum is

$$ \frac{\omega^h}{\omega} = \frac{1}{\omega h} \sqrt{\frac{20(2 - \cos(\omega h) - \cos^2(\omega h))}{16 + 13 \cos(\omega h) + \cos^2(\omega h)}} \quad (A.1) $$

The expansion $\cos(x) \sim 1 - x^2/2 + x^4/4! - x^6/6!$ is used, obtaining after some simple computations:

$$ \frac{\omega^h}{\omega} \sim \frac{1}{\omega h} \sqrt{\frac{30(\omega h)^2 - \frac{15}{2}(\omega h)^4 + \frac{11}{12}(\omega h)^6}{30 - \frac{15}{2}(\omega h)^2 + \frac{7}{8}(\omega h)^4}}, \quad (A.2) $$

which can be rewritten as

$$ \frac{\omega^h}{\omega} \sim \sqrt{\frac{N}{D}} \quad (A.3) $$

with N and D defined as follows:

$$ \begin{aligned} N &= 30 - \frac{15}{2}(\omega h)^2 + \frac{11}{12}(\omega h)^4, \\ D &= 30 - \frac{15}{2}(\omega h)^2 + \frac{7}{8}(\omega h)^4. \end{aligned} \quad (A.4) $$

Expression (A.3) can be written as

$$ \frac{\omega^h}{\omega} \sim \sqrt{\frac{1}{1 + \frac{D - N}{N}}} \quad (A.5) $$

and, using the expansions $\frac{1}{1+x} \sim 1-x$ and $\sqrt{1-x} \sim 1-x/2$, it gives rise to

$$ \frac{\omega^h}{\omega} \sim 1 + \frac{N - D}{2N}. \quad (A.6) $$

Finally, substituting the expressions (A.4) for N and D, we get

$$ \frac{\omega^h}{\omega} \sim 1 + \frac{(\omega h)^4}{1440}, \quad (A.7) $$

which reveals that the order of accuracy is 4. Fig. A.1 shows that, for low frequencies, the normalized discrete spectrum has the same behavior as the function $1 + (\omega h)^4/1440$ (recall: $\omega h = \pi n/n_{el}$).

<figure><img src="image_74.png" /><figcaption>Fig. A.1. Rod problem and consistent mass. Normalized discrete spectrum using quadratic NURBS compared with $1 + (\omega h)^4/1440$ for low frequencies.</figcaption></figure>
## A.2. Order of accuracy employing cubic NURBS and consistent mass

Using cubic NURBS, the normalized discrete spectrum is given by

$$ \frac{\omega^h}{\omega} = \frac{1}{\omega h} \sqrt{\frac{42(16 - 3 \cos(\omega h) - 12 \cos^2(\omega h) - \cos^3(\omega h))}{272 + 297 \cos(\omega h) + 60 \cos^2(\omega h) + \cos^3(\omega h)}} \quad (A.8) $$

Expanding $\cos(\omega h)$ and repeating the computations as before, we obtain:

$$ \frac{\omega^h}{\omega} \sim 1 + \frac{(\omega h)^6}{60,480} \quad (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

$$ \frac{\omega^h}{\omega} = \frac{1}{\omega h} \sqrt{\frac{2}{3}(2 - \cos(\omega h) - \cos^2(\omega h))}, \quad (A.10) $$

while with cubic NURBS we obtain:

$$ \frac{\omega^h}{\omega} = \frac{1}{\omega h} \sqrt{\frac{1}{15}(16 - 3 \cos(\omega h) - 12 \cos^2(\omega h) - \cos^3(\omega h))}. \quad (A.11) $$

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:

$$ \frac{\omega^h}{\omega} \sim 1 - \frac{(\omega h)^2}{8} \quad (A.12) $$

and for cubic NURBS

$$ \frac{\omega^h}{\omega} \sim 1 - \frac{(\omega h)^2}{6}. \quad (A.13) $$

As was already evident from Fig. 19, by increasing the order $p$, a better order of accuracy is not achieved. For “row sum” lumped mass, it is always equal to 2. Finally, Figs. A.4 and A.5 confirm the validity of expressions (A.12) and (A.13), respectively.

<figure><img src="image_75.png" /><figcaption>Fig. A.2. Rod problem and consistent mass. Normalized discrete spectrum using cubic NURBS compared with $1 + (\omega h)^6/60,480$ for low frequencies.</figcaption></figure>
<figure><img src="image_76.png" /><figcaption>Fig. A.3. Rod problem and “row sum” lumped mass. Analytical versus numerical discrete spectra computed using quadratic and cubic NURBS.</figcaption></figure>

<figure><img src="image_77.png" /><figcaption>Fig. A.4. Rod problem and “row sum” lumped mass. Normalized discrete spectrum using quadratic NURBS compared with $1 - (\omega h)^2/8$ for low frequencies.</figcaption></figure>

<figure><img src="image_78.png" /><figcaption>Fig. A.5. Rod problem and “row sum” lumped mass. Normalized discrete spectrum using cubic NURBS compared with $1 - (\omega h)^2/6$ for low frequencies.</figcaption></figure>
## References

[1] J.H. Argyris, The computer shapes the theory. Paper read to the Royal Aeronautical Society on May 18, 1965 (Lecture Summary: Aeronaut. J. Royal Soc. 69 (1965) XXXII).

[2] J.H. Argyris, Energy theorems and structural analysis. Part I. General theory, Aircraft Engrgr. 26 (1954) 347–356, 383–387, 394; 27 (1955) 42–58, 80–94, 125–134, 145–158; Part II (with S. Kelsey): Applications to thermal stress analysis and to upper and lower limits of St. Venant torsion constant, Aircraft Engrgr. 26 (1954) 410–422.

[3] J.H. Argyris, S. Kelsey, Energy Theorems and Structural Analysis, Butterworths, London, 1960.

[4] J.H. Argyris, S. Kelsey, Modern Fuselage Analysis and the Elastic Aircraft, Butterworths, London, 1963.

[5] ARPACK. Available from: <http://www.caam.rice.edu/software/arpack/>.

[6] J.K. Benninghof, R.B. Lehoucq, An automated multilevel substructuring method for eigenspace computations in linear elastodynamics, SIAM J. Sci. Comput. 25 (2004) 2084–2106.

[7] T.D. Blacker. CUBIT Mesh Generation Environment Users Manual, vol. 1, Technical Report, Sandia National Laboratories, Albuquerque, NM, 1994.

[8] L. Brillouin, Wave Propagation in Periodic Structures, Dover Publications, Inc., 1953.

[9] R.D. Buehrle, G.A. Fleming, R.S. Pappa, F.W. Grosveld, Finite element model development for aircraft fuselage structures, in: Proceedings of XVIII International Modal Analysis Conference, San Antonio, TX, 2000.

[10] D. Bushnell, Stress, stability and vibration of complex, branched shells of revolution, Comput. Struct. 4 (1974) 399–435.

[11] D. Bushnell, Computerized Buckling Analysis of Shells, Martinus Nijhoff Publishers, Dordrecht, 1985.

[12] A.K. Chopra, Dynamics of Structures. Theory and Applications to Earthquake Engineering, second ed., Prentice-Hall, Upper Saddle River, New Jersey, 2001.

[13] J. Chung, G.M. Hulbert, A time integration algorithm for structural dynamics with improved numerical dissipation: the generalized-α method, J. Appl. Mech. 60 (1993) 371–375.

[14] F. Cirak, M. Ortiz, Fully $C^4$-conforming subdivision elements for finite deformation thin shell analysis, Int. J. Numer. Methods Engrg. 51 (2001) 813–833.

[15] F. Cirak, M. Ortiz, P. Schröder, Subdivision surfaces: a new paradigm for thin shell analysis, Int. J. Numer. Methods Engrg. 47 (2000) 2039–2072.

[16] F. Cirak, M.J. Scott, E.K. Antonsson, M. Ortiz, P. Schröder, Integrated modeling, finite-element analysis, and engineering design for thin-shell structures using subdivision, Computer-Aided Des. 34 (2002) 137–148.

[17] R.W. Clough, J. Penzien, Dynamics of Structures, McGraw-Hill, New York, 1993.

[18] L. Couchman, S. Dey, T. Barzow, ATC Eigen-Analysis, STARS/ARPACK/NRL Solver, Naval Research Laboratory, Englewood Cliffs, NJ, 2003.

[19] G. Engel, K. Garikipati, T.J.R. Hughes, M.G. Larson, L. Mazzei, Continuous /discontinuous finite element approximations of fourth-order elliptic problems in structural and continuum mechanics with applications to thin beams and plates, and strain gradient elasticity, Comput. Methods Appl. Mech. Engrg. 191 (2002) 3669–3750.

[20] G.E. Farin, NURBS Curves and Surfaces: From Projective Geometry to Practical Use, A.K. Peters, Ltd., Natick, MA, 1995.

[21] I. Fried, D.S. Malkus, Finite element mass matrix lumping by numerical integration with no convergence rate loss, Int. J. Solids Struct. 11 (1975) 461–466.

[22] R.H. Gallagher, A Correlation Study of Methods of Matrix Structural Analysis, Pergamon, Oxford, 1964.

[23] F.W. Grosveld, J.I. Pritchard, R.D. Buehrle, R.S. Pappa, Finite element modeling of the NASA Langley Aluminum Testbed Cylinder. in: 8th AIAA/CEAS Aeroacoustics Conference, Breckenridge, CO, 2002, AIAA 2002-2418.

[24] H. Hilber, T.J.R. Hughes, Collocation, dissipation, and overshoot for time integration schemes in structural dynamics, Earthquake Engrg. Struct. Dyn. 6 (1978) 99–117.

[25] H. Hilber, T.J.R. Hughes, R. Taylor, Improved numerical dissipation for time integration algorithms in structural dynamics, Earthquake Engrg. Struct. Dyn. 5 (1977) 283–292.

[26] T.J.R. Hughes, The Finite Element Method: Linear Static and Dynamic Finite Element Analysis, Dover Publications, Mineola, NY, 2000.

[27] T.J.R. Hughes, J.A. Cottrell, Y. Bazilev, Isogeometric analysis: CAD, finite elements, NURBS, exact geometry, and mesh refinement, Comput. Methods Appl. Mech. Engrg. 194 (2005) 4135–4195.

[28] L. Meirovitch, Analytical Methods in Vibrations, The MacMillan Company, New York, 1967.

[29] I. Miranda, R.M. Ferencz, T.J.R. Hughes, An improved implicit–explicit time integration method for structural dynamics, Earthquake Engrg. Struct. Dyn. 18 (1989) 643–653.

[30] D.C. Montgomery, G.C. Runger, N.F. Hubele, Engineering Statistics, third ed., Wiley, New York, NY, 2003.

[31] E. Onate, M. Cervera, Derivation of thin plate bending elements with one degree of freedom per node: a simple three node triangle, Engrg. Comput. 10 (1993) 543–561.

[32] E. Onate, F. Zarate, Rotation-free triangular plate and shell elements, Int. J. Numer. Methods Engrg. 47 (2000) 557–603.

[33] R. Phaal, C.R. Calladine, A simple class of finite-elements for plate and shell problems: 1. Elements for beams and thin flat plates, Int. J. Numer. Methods Engrg. 35 (1992) 955–977.

[34] R. Phaal, C.R. Calladine, A simple class of finite-elements for plate and shell problems: 2. An element for thin shells with only translational degrees of freedom, Int. J. Numer. Methods Engrg. 35 (1992) 979–996.

[35] L. Piegl, W. Tiller, The NURBS Book (Monographs in Visual Communication), second ed., Springer-Verlag, New York, 1997.

[36] K.S. Pister, R.R. Reynolds, K. Willam (Eds.), Proceedings of Fenomech'78 – 1st International Conference on Finite Elements in Nonlinear Mechanics, University of Stuttgart, August 30–September 1, 1978. Computer Methods in Applied Mechanics and Engineering, vol. 17–18, 1979, 720 pp.

[37] S. Prudhomme, F. Pascal, J.T. Oden, A. Romkes, A priori error estimate for the Baumann–Oden version of the discontinuous Galerkin method, Comptes Rendus de l'Académie des Sciences I, Numer. Anal. 332 (2001) 851–856.

[38] D.F. Rogers, An Introduction to NURBS With Historical Perspective, Academic Press, San Diego, CA, 2001.

[39] G. Strang, G.J. Fix, An Analysis of the Finite Element Method, Prentice-Hall, Englewood Cliffs, NJ, 1973.

[40] C. Truesdell, The computer: ruin of science and threat to mankind, in: An Idiot's Fugitive Essays on Science: Methods, Criticism, Training, Circumstances, Springer-Verlag, New York, 1984, pp. 594–631.