NURBS-Based Collocation Methods for the Structural Analysis of Shells of Revolution

In this work we present a collocation method for the structural analysis of shells of revolution based on Non-Uniform Rational B-Spline (NURBS) interpolation. The method is based on the strong formulation of the equilibrium equations according to Reissner-Mindlin theory, with Fourier series expansion of dependent variables, which makes the problem 1D. Several numerical tests validate convergence, accuracy, and robustness of the proposed methodology, and its feasibility as a tool for the analysis and design of complex shell structures.


Introduction
The application of differential quadrature or collocation methods for the solution of engineering boundary value problems dates back to the initial works by Bellmann and Casti [1,2], who are credited to have introduced the term differential quadrature operator, as the counterpart of quadrature operator for the approximation of a differential operator, i.e., the strong form of the equation; and having applied this strategy to first-order ordinary non-linear differential equations.The terminology of collocation method is, to the best of the authors' knowledge, associated with the work by Petersen [3], and subsequently with the work by Karpilovskaya [4], Vainikko [5,6] and Kaspshitskaya [7].The two methods share substantial similarities, namely the main idea of solving the strong form of the problem with a proper choice of discrete function spaces of approximation.By introducing a proper interpolation of the dependent variables, which is substituted in the governing differential equations, a linear set of equations for the unknown variables can be written by enforcing (i.e., collocating) the discretized form of the governing equations at a set of locations on the domain, known as collocation points.
The method appears to have gained more interest in the scientific research community in the last two decades, rather than immediately following the presentation in the scientific literature.Nowadays, there is discussion as to whether this methodology may become a valid and general purpose alternative to more established methods such as the finite element method [8].Without trying to give a general answer to this question, and without trying to provide a detailed literary review of the many applications of the collocation or differential quadrature method in past works, the aim of the present communication is to present an innovative formulation of the collocation method for the analysis of shells of revolution based on Non Uniform Rational B-Spline (NURBS) functions.
The collocation method is ideal in application to structural theories and has found several instances in the scientific literature [9].With reference to the analysis of shell structural elements, Mirfakhraei and Redekop [10] applied a Lagrangian differential quadrature method to the study of buckling phenomena of orthotropic thin shells of revolution; Wu et al [11,12] presented solutions to problems of axisymmetric bending of shells of revolution, introducing the so-called generalized differential quadrature rule.Viola and Artioli [13,14] and Artioli et al. [15,16] proposed a Lagrange differential quadrature solution for complex rotational shells with Fourier series expansion of dependent variables, both for the free vibration analysis, and for the static analysis of shells of revolution subject to unsymmetrical loading.The application of a collocation or differential quadrature method using Spline interpolants had been already proposed by Bellman and co-workers to solve general fourth-order differential equation problems [17], and later applied by Zhong and Guo [18,19] to static and buckling problems for Kirchhoff plates and to the non-linear free vibration analysis of straight Euler beams, respectively.Dade Huang proposed a Spline based collocation method for laminated plates and shells [20], while Viswanathan introduced the collocation method with Spline basis function in application to shear deformable laminated plates and shells [21][22][23].
The recent exponential growth of the isogeometric analysis [24][25][26][27][28][29] has raised attention on collocation methods based on NURBS and/or T-Splines trial functions.Given their higher-order smoothness and their favorable approximation properties, NURBS basis functions constitute an efficient interpolation tool for collocating differential operators, i.e., for the approximation of the strong formulation of PDEs; furthermore they apply even to irregular or complex geometric domains.The newly developed IGA Collocation methods [8,[30][31][32][33], set in the realm of weighted residuals [35], have been applied to a variety of solid, structural mechanics and engineering problems, with excellent performances in terms of efficiency, accuracy and robustness if compared to standard Galerkin methods.
The objectives of this work are the development and validation of the proposed methodology as an efficient, accurate numerical tool for the analysis and design of shells of revolution.The paper is organized as follows.In section 2 we formulate the differential problem of free vibration and of static equilibrium for a shell of revolution in the framework of Reissner-Mindlin theory.In this section, consistent cross-reference to Appendix A-Fundamental relationships for shells of revolution-and Appendix B-Equilibrium and boundary conditions operators-is made to make treatment more concise.In Section 3 some introductory elements on B-Splines are recalled from classic geometrical interpolation theory.Section 4 presents the basic one-dimensional NURBS interpolation basis functions as a result of certain projective transformations of B-Spline function spaces and its application to the NURBS-based collocation method, object of the present work.Section 5 collects an extensive set of numerical tests focused on the analysis of the convergence properties of the proposed method, as well as the validation in terms of accuracy and robustness by comparison with existing results.Last, the collocation method is applied to a problem of design for a complex shell structure, showing its feasibility as a numerical tool for engineering applications.Conclusions and future work perspectives are drawn in Section 6.

Geometry
We consider a shell of revolution defined by its mid-surface, obtained by rotating a plane generator curve around an axis to form a closed surface, as shown in Figure 1a.Curvilinear coordinates on the mid-surface are chosen to be associated to two families of mutually orthogonal curvilinear lines over the mid-surface, i.e., meridians and parallels, respectively.In particular, z is the coordinate measured along the axis of rotation, any value of z indicates a fixed parallel, whereas θ is the circumferential angle, indicating a given meridian.The third orthogonal coordinate is ζ, the distance of a point from the mid-surface, measured along n, the outward unit normal to the mid-surface.The shell thickness will be denoted by t.Thus, the mid-surface can be described by the bounded open set: and its boundary For subsequent developments, it is here assumed that description of shell geometry is given in terms of the parallel radius equation R = R(z), which defines the shape of the meridian curve, plus the thickness variation t = t(z), for z ∈ [z i , z f ].Nonetheless, it is pointed out that this hypothesis is not mandatory: R(z) and t(z) need not be specified analytically; i.e., they can be obtained from interpolated data, for instance with the use of NURBS basis functions [36].According to the above geometry description, metrics of the shell of revolution mid-surface are resumed in Appendix A.

Displacement Field
Following the classical Reissner-Mindlin theory [37,38], and referring to Figure 1a, the displacements of a material point of the shell are u, v, and w, where u is the translation along t 1 , the unit vector tangent to the mid-surface and directed along the local meridian, v is the translation along t 2 , the unit vector tangent to the mid-surface and directed along the local parallel, and w is the displacement along n, orthogonal to the mid-surface.Moreover, β 1 and β 2 represent the rotation of n about meridian t 1 and circumferential t 2 direction, respectively.In accordance with the aforementioned theory, the components u and v vary linearly with respect to ζ, while w, β 1 , and β 2 are constant through the thickness.The degrees of freedom are collected in the algebraic vector: depending on the mid-surface coordinates z and θ.

Internal Actions
In the present theory, equilibrium equations are written in terms of internal actions, or stress resultants and couples per curvilinear coordinate line unit length, defined as the integrals of the stress components through the shell thickness [37].In Figure 1b, stress resultants and couples are represented acting on a mid-surface element, identified by two couples of curvilinear coordinate lines, having in-plane outward normal unit vectors t 1 and t 2 , respectively.These internal actions can be grouped as membrane normal N 11 , N 22 , [resp.shearing N 12 ] forces; flexural bending M 11 , M 22 , [resp.twisting M 12 ] couples; and transverse shear forces Q 11 , and Q 22 , respectively.The set of stress resultants and couples are collected in the algebraic vector: Finally, q(z, θ) = [q u q v q w 0 0] T are surface distributed load components per unit mid-surface area corresponding to the degrees of freedom u.Surface distributed bending loads are not considered in the present development.

Equilibrium
Following classical arguments [38], the free vibration problem, as well as the static equilibrium problem for a shell of revolution, are formulated as a set of five differential equations for the internal actions, a format which is classified as statically indetermined [37].Nevertheless, the above problem may be formulated as five differential equations for displacements and rotations plus boundary conditions, introducing proper shell-like strain measures, defined through a strain-displacement linear differential operator and choosing an isotropic linear elastic constitutive law between internal actions and strain components (see Appendix A).The strain components are conjugate to internal actions S, and are collected in the algebraic vector: Applying Fourier series expansion of dependent variables (3), (4), and (5) with respect to the circumferential angle θ, the latter form can be reduced to a series of sets of five linear differential equations for the displacement harmonic amplitudes, on a one dimensional domain: the meridian curve.Since this treatment is quite consolidated [37][38][39][40], and the present work focuses on the Fourier series form of the equilibrium operator, for conciseness, we hereby present the displacement formulation in the so called harmonic form, postponing notations and the above recalled intermediate derivations to Appendices A and B.

Free Vibration Problem
The free vibration problem is governed by the classical linear eigenvalue problem: Find (u (j) , ω (j) ) s.t. for any natural number j, L (j) (u (j) ) = (ω (j) ) 2 Mu (j) , ∀z ∈ (z i , z f ), where (u (j) , ω (j) ) is the couple formed by the jth eigen-mode harmonic amplitude (in the sense specified in ( 35)), and the associated natural frequency.The inertia operator M is here taken as: where ρ is the material mass density [38].

Static Equilibrium Problem
The static equilibrium problem is governed by the following linear elliptic problem: The equilibrium operator L (j) and the boundary conditions operator B (j) appearing in ( 6) and ( 7) are reported in Appendix A, depending on the equation R = R(z), i.e., on the meridian shape, and on the relative metric quantities.

B-Spline Basis Functions
Non-Uniform Rational B-Splines -NURBS -are obtained from linear combination of B-Splines basis functions locally defined over patches of intervals [41].A B-Spline of degree p is defined through an ordered set of non-decreasing coordinates in the parametric space: the knot vector Ξ = {ξ 1 = 0, ξ 2 , . . ., ξ n+p+1 = 1}, where p represents the polynomial degree of the interpolation, and n is the number of basis functions associated to the knot vector.Consequently, it results that p = 0, 1, 2, 3 . . .refers to constant, linear, quadratic, cubic,... piecewise polynomials, respectively.An open knot vector has the first and the last knots repeated p + 1 times ensuring that the basis functions are interpolatory at the extrema of the parametric space [41].In the following, we assume uniform, i.e., equally spaced, open knot vectors.
The B-Spline basis functions N i,p can be defined starting with piecewise constants (p = 0) by means of a recursive formula For p = 1, 2, 3, ... it results:

, otherwise
Accordingly, basis functions of order p have (p − 1) continuous derivatives.If a knot is repeated k times, then the number of continuous derivatives decreases by k.A knot with multiplicity p, implies the basis functions to be interpolatory at that knot.

One Dimensional NURBS Interpolation
A non-uniform 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 [41,42].NURBS can be used to construct a vast gallery of geometrical entities, starting, for instance, from a discrete set of measured data [36].The projective transformation of a B-Spline curve yields a rational polynomial curve.Inductively, the order p of a NURBS basis is identified as the order of the polynomial basis from which it was generated by projection.
A NURBS interpolation in R d is constructed starting from a data set B w i (i = 1, ..., n) or control points, for a B-Spline curve in R d+1 with knot vector Ξ.The control points for the NURBS curve are given by [41]: where (B i ) k is the kth 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: Their first and second derivatives are given, respectively, by: and A pth degree NURBS curve is defined by:

NURBS-Based Collocation Method
The main idea in the collocation method is to approximate the strong form of the equilibrium equations, without addressing a weak formulation typical of the finite element approach [8].With reference to problem (6) or (7), the method consists in the approximation of the solution u (j) , j = 0, 1, 2, ..., as: where the interpolation functions R i,p (ξ) follow (10).The discrete form of ( 6)-( 7) is thus obtained substituting ( 14) into the equilibrium operator L (j) and into the boundary conditions operator B (j) , respectively.This leads to a set of 5n linear equations for the control values corresponding to (14).These unknown control values are obtained by enforcing the resulting discrete governing equations at n prescribed points in the parametric space, ξ c (c = 1, ..., n), defined collocation points: • Static equilibrium problem Solution is thus obtained by solving for the unknown control values ũ(j) i .Collocation points are a set of natural coordinate values ξ c , c = 1, ..., n, related to the knot vector.In the present case we will adopt the so called Greville nodes [36,43], i.e.,: The choice for collocation points is generally guided by an accuracy/efficiency criterion.The idea is to have a set of points for which known solutions of certain differential operators are exactly recovered by the numerical scheme [17,[44][45][46].The choice for collocating governing equations as in Equation ( 17) has proved efficient in previous NURBS-based collocation approaches [8,32,33], as it grants high-order accuracy and smoothness of the parametrization of the unknown with excellent convergence properties.Noteworthy, according to Johnson [46], Greville collocation points are associated to the property of linear precision [43], in the sense that if these are used to compute a B-Spline approximation of a straight line of the form l(u) = au + b, the corresponding B-Spline curve is a straight line.
Given the interesting performance of standard Lagrange interpolation spaces, for the implementation of 1D collocation methods for shells of revolution, reported by previous researches [13][14][15][16], in what follows, we explore the superior accuracy of the NURBS approximation in conjunction with the collocation approach.The attractiveness of this approximation is twofold: (i) the NURBS high regularity suggests the possibility of obtaining smooth solutions even with low order approximation and low number of degrees of freedom, (ii) the capability of describing complex geometries is interesting especially from a structural design point of view, when constructing geometrical models from measured data is required.
The discretization and solution of problem Equations ( 6) or ( 7) is programmed as follows: 1. Define NURBS order p, number of control points n, build open knot vector Ξ, build NURBS basis and derivatives; 2. Compute collocation points ξ c , c = 1, 2, ..., n; 3. Parametrize curvilinear space coordinate z = [z i , z f ] as z = z(ξ), and relevant geometric quantities, i.e., R(z(ξ)), A 1 (z(ξ)), sin φ(z(ξ)); 4. Collocate above mentioned quantities at collocation points; 5. Approximate equilibrium, mass matrix, and boundary conditions operators: L (j) (z(ξ)), M(z(ξ)), and B (j) (z(ξ)), respectively, and applied load distribution q (j) (z(ξ)), through NURBS interpolants; As it is well known, such a discretization procedure does not lead to symmetric banded linear sets of equations as it is the case with finite element procedures.In particular, the coefficient matrix is sparse and nonsymmetric [44]; accordingly, factorization methods are generally not the most efficient, especially when dealing with large scale simulations.It is therefore recommended to resort to iterative methods such as conjugate gradient methods and residual based methods; is nonetheless noted that, albeit nonsymmetric, the pattern of coefficient matrices in displacement-based collocation approaches results competitive with respect to Galerkin finite element approaches which as known lead to symmetric systems [8].

Numerical Results
The present section is dedicated to numerical tests validating the NURBS collocation method for shells of revolution presented in the previous section.In subsection 5.1, we first present the asymptotic analysis of the fundamental eigenfrequency for shells of revolution of different Gaussian curvature (see (34)), an evidential benchmark test for assessing robustness with respect to membrane and shear locking in rotational shell numerical formulations.Secondly, we address the classical free vibration problem for different shells of revolution and different boundary conditions, comparing results with reference solutions.In subsection 5.2, a set of static problems are taken into consideration, in order to address the question of accuracy and convergence with respect to the order of interpolation and discretizaion.Finally, a simple technical problem is considered, validating the applicability of the present approach to problems of shell analysis and design, for shells with complex geometry.

Asymptotic Analysis of Fundamental Frequency
Following [47][48][49], in this section we study the dependence of the lowest (fundamental) frequency on the thickness t as this quantity approaches zero.Three representative cases of Parabolic (zero Gaussian curvature), Elliptic (positive Gaussian curvature), and Hyperbolic (negative Gaussian curvature) shells of revolution are considered.The shells have clamped boundary conditions, in order to study the inhibited pure bending behavior [49].The three cases reported in the following are characterized, respectively, by meridian equation of: • Parabolic shape: • Elliptic shape: with a 2 = 25 9 H 2 , b 2 = H 2 ; • Hyperbolic shape: with and H = 1.Material constants for the numerical simulations are E = 2.1 × 10 11 Pa, ν = 0.3, ρ = 7868 Kg/m 3 .The scope of this test is to check the asymptotic behaviors of the square of the fundamental, i.e., lowest, eigenfrequency ω 2 0 when computed with the NURBS collocation method.In passing, it is recalled that the lowest eigenfrequency ω 2 0 does not necessarily correspond to harmonic number j = 0, nor any specific a-priori-known value j.Expected asymptotic rates are [47,50]: • Elliptic shell: • Hyperbolic shell: In Figure 2 we report the logarithmic dependence of ω 2 t vs. thickness ratio t/H, for the three considered geometries.Accuracy in the low thickness range requires a fine discretization and high NURBS interpolation order, namely n = 101 and p = 11.In Table 1 we report the eigenvalues at specified thickness values, delimiting subsequent thickness branches in the logarithmic scale.The slopes of these branches are in agreement with the theoretical prediction on asymptotic behavior of the fundamental frequency Equations ( 21)- (23).The above results prove the present method is capable of efficient treatment of situations in which geometry, boundary conditions and regularity of applied loads induce a so called intermediate state for which nor the bending strain energy or the membrane strain energy dominates on one another in the asymptotic regime, and interior or boundary layers concentrate the majority of the strain energy.Such situations are prone to numerical locking and strong instabilities in the solution [34,48].In this regard, a necessary condition for a numerical scheme to be locking-free is given by Equations ( 21)- (23), which are satisfied by the present method.

Clamped Hemispherical Shell in Free Vibration
A clamped hemispherical dome with radius a and uniform thickness t (a/t = 100), see Figure 3, is considered [51,52].The values of the non-dimensional frequency parameter λ = aω (j) (ρ/E) 1/2 are reported in Table 2, for the first six modes and for different values of the harmonic number j. Reference results obtained with a Lagrange collocation method [14], and those from a Spline finite element solution by Luah and Fan [53] are reported for validation.The present frequency parameter values are obtained selecting n = 35, p = 4 for the lowest mode frequencies, up to n = 85, p = 8 for higher modes.Noteworthy, the closed apex condition required by this type of shell [16,37] can be modeled by the present approach straightforwardly, as it implies enforcement of special homogeneous Dirichlet boundary condition at z = 0, depending on the harmonic number j: u (1) + v (1) = 0, w (1) = 0, β (1) The present results are in agreement with the reference results, confirming the accuracy of the proposed technique, even for cases in which shear locking is expectable [54].The method proves computationally efficient as even in the higher mode regime, the size of the eigenvalue problem is smaller than that associated to the reference results.

Simply supported circular plate in free vibration
A circular plate of radius a and thickness t, simply supported along the boundary, is considered.For conciseness, we shall refer to the shell geometry and boundary conditions as shown in Figure 7 even if related to the numerical test reported in Section 5.2.2.The non-dimensional parameter λ = [12(1 − ν 2 )ρa 4 (ω (j) ) 2 /(Et 2 )] 1/2 computed with the present method is reported in Table 2, in agreement with reference values obtained in [14], and with those available from a Spline finite element approach by Luah and Fan [53].The results confirm the interesting capabilities of the proposed approach for the free vibration analysis of rotational shells with different shapes and boundary conditions at a low computational cost.Numerical simulations show monotonic convergence to the computed eigenfrequencies, for increasing interpolation order p up to p = 5, and discretization order n up to n = 35.No spurious modes arise.Table 2. Clamped hemispherical shell in free vibration.Non-dimensional frequency parameter λ = aω (j) (ρ/E) 1/2 .Reference results from [53] (parentheses) and from [14] (brackets).

Harmonic Number j
Mode

Cantilever Cylindrical Shell under Tip Load
A circular cylindrical shell, clamped at one end, with a uniform radial line load Q = 1 kN/m acting on the free edge, is considered, see Figure 4.The geometric data are R = 5 m, t = 0.01 m, H = 6 m, the Young's modulus and the Poisson ratio are, respectively, E = 10 7 kN/m 2 , and ν = 0.3.Axial symmetry permits to address the equilibrium problem for the j = 0 harmonic only.The line load is taken into account by modifying the free boundary conditions at z = 0, Equations ( 59)-( 61), to guarantee that Q 11 equals the applied force on that boundary.
The present classical benchmark is here adopted as a tool for calibrating the model with respect to parameters n and p, i.e., mesh discretization (h−refinement) and order of interpolation (p−refinement), and for a given accuracy level.To this end, a reference solution proposed by Timoshenko and Woinowsky-Krieger [55] for semi-infinite circular cylindrical shells in axisymmetrical bending, is here taken into consideration for comparison and accuracy assessment.The structural behavior, which shifts from shear-flexural near the loaded edge to membrane dominated at the clamped boundary, is correctly reproduced by the NURBS-based collocation approach, as can be seen in Figure 5a,b, where M 11 and Q 11 are plotted, respectively, for an overkilling discretization with n = 250 and p = 7 with a relative error norm below 10 −12 , together with values computed through the reference solution at selected locations.The results are in agreement, plus the static boundary condition is exactly satisfied in solution.Moreover, in order to assess the method convergence properties, we have plotted the H 1 -error norms, of the tip transverse displacement w(0), Figure 6a, and of bending moment M 11 (0), Figure 6b, for odd and even NURBS order p, and for increasing discretization refinements n.The asymptotic tendency of the error norm is consistent with analogous results reported in [31], where the static problem for Timoshenko spatial beams is addressed with a NURBS based collocation method.In particular, the logarithmic error is of the order n −p for any duplets of orders (p, p + 1), with even p, as theoretically predicted [32].

Simply Supported Circular Plate Subject to Uniform Pressure
We consider a simply supported circular plate under uniform normal load q w , Figure 7.This benchmark has a closed-form solution, cfr.[56].The material parameters are E = 10 6 kN/m 2 , ν = 0.0, the plate radius is a = 1 m, the thickness is t = 10 −3 m, and the applied surface load is q w = 1 kN/m 2 .The load presents axial symmetry and requires to solve (7) for the j = 0 harmonic only; closed apex conditions plus straight meridian geometry is fully taken into account by the present approach as the former amounts to special homogeneous Dirichlet boundary conditions (24), the latter is contemplated by ( 30)- (31).We report in Table 4 the results obtained with the present approach, compared with an analytical solution [37], and with a hybrid mixed rotational finite element solution [57].A space discretization characterized by n = 7 and p = 4 is sufficient to obtain accurate results.Resorting to the NURBS basis structure, post-processed quantities such as the bending and twisting moments in Figure 8a, and the rotation degrees of freedom in Figure 8b, can be plotted.

Table 3.
Simply supported circular plate in free vibration.Frequency parameter λ = [12(1 − ν 2 )ρa 4 (ω (j) ) 2 /(Et 2 )] 1/2 .Reference results from [14] (parentheses) and from [53] (brackets).This numerical test presents a spherical shell of radius a = 1 and uniform thickness t = 0.01, see Figure 9.The shell is loaded by a uniform internal pressure q w = 1.The Young modulus is E = 2.1 × 10 11 and the Poisson ratio is ν = 0.3.Due to symmetry, only the upper half shell is modeled, with closed apex conditions at z = 0 and sliding edge boundary condition at z = R.The solution for displacement and strain components is expected to in a pure membrane dominated regime with zero shear/bending deformation.Henceforth, given the aim of the present test, no reference solution is here taken into account.The NURBS collocation is applied with n = 35 and p = 5, sufficient to grant the exact solution reported in the following.The results shown in Figure 10a,b for displacement and strain components are in agreement with theory, in particular zero tangential displacements v as well as null membrane shearing strains ε 12 are correctly recovered, for a pure membrane equilibrium state of the spherical shell characterized by ε 11 = ε 22 .It is readily checked that curvature variations and shear strains are identically zero as expected.

Hyperboloid of Revolution under Wind Load
The example presents a doubly curved hyperboloid of revolution, a typical cooling tower structure, Figure 11, subject to wind pressure.Geometrical data for the shell are, a = 84 ft, T = 60 ft, H = 330 ft, t = 7 in.The Young modulus is E = 432 × 10 6 lb/ft 2 , Poisson ratio is ν = 0.15.The surface load q w is constant along the tower height, and is not symmetric with respect to the axis of revolution; in particular, q w varies with the circumferential coordinate θ as reported in Figure 12 [57].The distributed load is expanded in Fourier series taking into account the first ten harmonic components j = 0, 1, ..., 9.The present analysis aims at determining both displacements and stress resultants for the structure, which are respectively plotted in Figures 13 and 14, where a comparison with different reference solutions available in the literature is carried out [58], together with a reference solution computed with SHORE III, a ring finite element code developed by [59].The proposed methodology yields reliable results for n = 35 and p = 8 even for the case of a doubly curved rotational shell with unsymmetrical loading, which produces complex stress resultants and couples patterns.The numerical test confirms the present method capability of analyzing real shell structures under complex loading conditions with favorable accuracy.blue solid line-present approach vs. blue bullets-reference results [58]; membrane stress resultant N 22 : red solid line -present approach vs. red bullets-reference results [58].(b) Bending moment M 11 : blue solid line-present approach vs. blue bullets-reference results [59]; bending moment M 22 : red solid line -present approach vs. red bullets-reference results [59].

Torispherical tank subject to uniform internal pressure
The example presents a technical problem for a so called compound shell, namely a shell of revolution obtained by the smooth blending of different shell segments, each having a particular meridian curve geometry [37], see Figure 15.The shell comprises a closed-apex spherical cap OA 1 of radius a 1 = 172.8981in.and length H 1 = 18.675 in., a torus segment A 1 A 2 of radius a 2 = 32.74in.and length H 2 = 29.2in., and a cylindrical portion with radius a 3 = 96.1 in.and length H 3 = 24 in.; the thickness t is uniformly equal to 0.20 in.The Young modulus is E = 3.0 × 10 4 psi, the Poisson ratio is ν = 0.3.In order to model the pre-buckling elastic stress pattern, which is known to induce particularly severe compressive circumferential stresses in the toroidal region, the shell is loaded by a uniform internal pressure q w = 150 psi [37].The geometry in the present analysis is constructed using three different NURBS control polygons, each one pertaining to a single shell segment, and each set with an equal number of basis functions.The three subdomains are then patched together, granting regularity of the meridian at the junctions and C 0 continuity of displacements and rotations, a procedure referred to as domain decomposition [44,45].Nonetheless, it is pointed out that it is possible to use a single basis to describe the entire geometry and to collocate governing equations over the single control polygon.NURBS basis for collocation is characterized by n = 55 on the spherical and toroidal segments, and n = 7 on the cylindrical region.The order of the interpolants is p = 5.The reference solution for this benchmark is provided by the finite element code SHORE III, adopting a fine mesh discretization with 128 quartic displacement based elements on both spherical and toroidal parts, and with 8 elements on the cylindrical segment.Results and comparison are plotted in Figure 16, reporting displacements u and w, and in Figure 17 for rotation β 1 .Figure 18a [resp.18b] show the axial σ 11 and circumferential σ 22 normal stress components at inner (ζ = −t/2) and outer (ζ = t/2) surface of the shell [60]: obtained by post processing of stress resultants and couples.Remarkably good agreement is observed between present and reference results, confirming the capability of the proposed technique to study complex technically relevant problems for shells of revolution, as the aforementioned pre-buckling analysis of a cylindrical tank with torispherical head. .Torispherical tank subject to uniform internal pressure.Displacement components: u blue solid line -present approach vs. blue bullets -reference solution [53], w red solid line -present approach vs. red bullets -reference solution [59].red solid line-present approach vs. red bullets-reference solution [59].(b) Circumferential normal stress σ 22 -inner (ζ = −t/2) blue solid line-present approach vs. blue bullets-reference solution [59], outer (ζ = t/2) red solid line -present approach vs. red bullets-reference solution [59].

Conclusions
In this paper we have presented a NURBS-based collocation method for the elastic analysis of shells of revolution.Both the free vibration and the static equilibrium problems have been investigated in the realm of Reissner-Mindlin theory.Use of Fourier series expansion of dependent variables with respect to the circumferential angle permits to solve the problem as a series of one-dimensional differential equation sets.A NURBS-based collocation method has been presented exploiting the superior behavior of such interpolation methodology which grants higher order accuracy and regularity.The proposed method has been tested on different free vibration and static equilibrium benchmark problems, showing its versatility and accuracy in computing free vibration modes and frequencies, as well as stresses and displacement patterns.Comparison with existing results seems encouraging and confirm the interesting approximation properties of NURBS basis functions which grant avoidance of numerical problems such as locking and instability.The robustness shown by the proposed methodology suggests implementation into a general purpose code for the structural analysis and design of complex shell structures.

Figure 1 .
Figure 1.Geometry of shells of revolution.

Figure 2 .
Figure 2. Asymptotic analysis of fundamental frequency.Asymptotic pattern of lowest eigenfrequency ω 0 with respect to thickness ratio t/H approaching zero.Different shells of revolution: elliptic -red dashed line with stars, parabolic -blue solid line with circles, hyperbolic -black dashed-dot line with bullets.

Figure 3 .
Figure 3. Clamped hemispherical shell in free vibration.Cross section.

Figure 7 .
Figure 7. Simply supported circular plate subject to uniform pressure.Cross section.

Figure 15 .
Figure 15.Torispherical tank subject to uniform internal pressure.Cross section.

Figure 16
Figure16.Torispherical tank subject to uniform internal pressure.Displacement components: u blue solid line -present approach vs. blue bullets -reference solution[53], w red solid line -present approach vs. red bullets -reference solution[59].

Table 1 .
Asymptotic analysis of fundamental frequency ω 2 0 with respect to thickness ratio t/H for elliptic, parabolic, and hyperbolic shells of revolution.Slope values for different branch values in the asymptotic regime.

Table 4 .
[57]ly supported circular plate subject to uniform pressure.Comparison of present results for w, β 1 , and M 11 at different locations z, for cp = 7, p = 4. Theoretical results from[56], reference results from[57].