Prediction of natural frequencies of laminated curved panels using refined 2-D theories in the spectral collocation method

Abstract This paper presents a versatile and efficientmodeling and solution framework for free vibration analysis of composite laminated cylindrical and spherical panels modeled according to two-dimensional equivalent singlelayer and layerwise theories of variable order.Aunified formulation of the equations of motion is adopted which can be used for both thin and thick structures. The discretization procedure is based on the spectral collocation method and is presented in a compact matrix form which can be directly and easily implemented. The convergence and accuracy of the proposed approach is evaluated for panels having different boundary conditions, thickness and shallowness ratios, and lamination layups.


Introduction
Accurate yet efficient modeling and solution of shell structures can be considered a classical problem in structural mechanics [1][2][3][4]. Over the last 50 years, researchers proposed many different approaches to face this problem [5][6][7], urged by the increasing adoption of curved panels as structural elements in a large variety of engineering applications. For example, cylindrical and spherical panels are important components of built-up aircraft and spacecraft structures. Shell structures are also widely applied in automotive, marine and civil engineering disciplines. In the last three decades, the design of shell structures has known new exciting opportunities by the usage of composite materials [8,9]. Composite laminated shells offer higher stiffness/strength to weight ratios than most metallic constructions and today they are fairly common in advanced engineering applications such as aerospace systems.
In service, such curved laminated panels are typically subjected to various dynamic loads, and it is crucial that their dynamic response is well predicted from the design stage onwards in order to assure their integrity and stability. One of the most important purpose in any linear dynamic analysis is the accurate computation of the natural frequencies of the system [10,11]. Broadly speaking, the accuracy of free vibration analysis of laminated shells mainly relies on two aspects. The first issue is related to the assumptions and simplifications adopted in the mathematical model of the structural component. The second aspect involves the method selected to solve the governing equations of the problem. A brief discussion of both aspects follows.
Models of curved panels based upon the three dimensional (3-D) theory of elasticity can be considered as the most accurate, since no overly simplified assumptions are introduced in describing the kinematics of deformations [12][13][14][15]. As such, 3-D models are suitable for shells of any thickness ratio (defined as the ratio between the thickness of the panel to the shortest of the span lengths or radii of curvature) and any shallowness ratio (defined as the ratio of the shortest span length to one of the radii of curvature), ranging from thin and shallow to thick and deep shells vibrating at low to high frequency. However, a full 3-D dynamic analysis, especially for composite shells, is rather complicated and time consuming. Therefore, socalled shell theories aimed at reducing the problem from three to two dimensions have been introduced by employing appropriate assumptions on the displacement behavior in the thickness direction.
Many different displacement-based two-dimensional (2-D) theories for laminated curved panels are available and a review of the extensive literature on this topic is beyond the scope of the present work [9]. Generally speaking, we can distinguish between equivalent single-layer (ESL) and layerwise (LW) theories [16], including or not shear deformation, rotary inertia and thickness stretching factors. ESL approaches assume a proper overall kinematic field throughout the thickness of the laminated structure, whereas an independent displacement field is postulated for each layer in a LW framework and appropriate continuity conditions are imposed at each layer interface. Typically, both ESL and LW kinematic fields are expressed as complete polynomial series expansion of the thickness coordinate and the highest power of the assumed polynomial set is generally referred to as the order of the theory. Owing to their intrinsic simplifications, 2-D theories provide reliable models for a limited range of thickness ratios, frequency values and through-the-thickness variation of material properties. The accuracy usually degrades as the wavelength of the vibration mode is of the order of magnitude of the panel thickness and as the variation of mechanical properties through the thickness direction increases like the case of sandwich panels. This loss of accuracy can be successfully compensated by using theories of higher order and/or relying on a LW approach, with the price of increasing the complexity and computational cost of the resulting models. Frequently, the knowledge in advance of the correct shell theory balancing the accuracy and computational burden for the specific problem under investigation can be a hard task. A unified modeling framework capable of tailoring the order and typology of the shell theory without the need of a new modeling effort each time would be highly desirable.
The other important aspect related to the computation of natural frequencies of laminated curved panels is the method adopted to solve the equations governing the dynamic problem. Again, many different approaches are available in the literature [9]. Exact solutions satisfying both the differential equations and the boundary conditions are possible only for a limited set of shell geometries, boundary conditions and lamination sequences. In most practical situations, one must rely on approximate methods. The most common and traditional approaches include the Ritz method and the finite element method (FEM). For simple structures, the Ritz method shows better convergence and less computational need than FEM. However, the FEM overcomes the limitations of the Ritz method in dealing with complicated boundary conditions and complex shapes. More recently, new emerging meshless methods are increasingly applied to the analysis of shell problems with the aim of eliminating some difficulties existing in FEM such as mesh distortion and remesh-ing [17]. Among recent solution methods, the generalized differential quadrature (GDQ) is also attracting more and more interest due to its simplicity and versatility [18][19][20].
The scope of this work is to present an advanced modeling and solution framework for the free vibration analysis of both thin and thick, deep and shallow laminated cylindrical and spherical panels with different combinations of boundary conditions. The range of applicability of the resulting tool is similar to a fully 3-D analysis, with also the possibility of using economical low-order theories when more refined models are not required. Therefore, the balance between accuracy and computational savings can be tailored on the specific application under study. The modeling aspect exploits the power and versatility of the unified formulation proposed by Carrera [21], which provides a smart way of handling arbitrary refinements of classical shell theories. The discretization of the problem to obtain an approximate solution of the natural frequencies is performed by the spectral collocation method, also called pseudospectral method [22,23], which is known to exhibit high rate of convergence and accuracy for differential problems without singularities. In so doing, relatively light discretized models are obtained which can be profitably used in extensive optimization and parametric studies.
As a final remark, it is noted that recent papers with some similarities with the present approach are available in the literature. In particular, worth mentioning are the works of [20], [24] and [25]. The first reference contains a comprehensive study of free vibrations of eight different shell structures using the GDQ method. The work is highly insightful but it is limited to models based on ESL theories. The paper by [24] applies a radial basis functions collocation method to perform a static and free vibration analysis of cylindrical and spherical laminated panels. Although the approximation of the problem is different, the discretization procedure is quite similar to what presented here. However, the analysis has the limitation of considering only models based on a first-order layerwise theory. Finally, the work reported in [25] uses a hierarchical trigonometric Ritz method to compute the natural frequencies of laminated shells. The approach is similar to that provided by [26][27][28] and can rely on both ESL and LW models of variable order. However, the analysis in the paper by [25] considers only panels having fully simplysupported edges. This work aims at overcoming the limitations of the above references by considering a comprehensive set of both ESL and LW shell theories of different orders and studying laminated panels with various combinations of classical boundary conditions.

Modeling framework
Let's consider the cylindrical and spherical laminated panels in Figure 1, which are composed of N ℓ layers of homogeneous orthotropic material. Each layer k has thickness h k and is numbered sequentially from bottom (k = 1) to top (k = N ℓ ) of the panel. The total thickness of the panel is h = ∑︀ N ℓ k=1 h k . The undeformed middle surface Ω k of each layer is described by the two orthogonal curvilinear coordinates α and β. Let z k denote the rectilinear local thickness coordinate in the normal direction with respect to Ω k . The components of the displacement field of layer k are indicated as u k α , u k β and u k z in the α, β and z directions, respectively.
According to 3-D elasticity and considering curved panels with constant curvature, the in-plane strains ε k p = where (2) R k α and R k β are the curvature radii of the α and β coordinate curves, respectively, at the generic point of the middle surface Ω k of the layer, and It is noted that H k α = 1 for cylindrical panels since R k α = ∞ and H k α = H k β for spherical panels since R k α = R k β . Note also that when 1/R k α = 1/R k β = 0, the above relations degenerate to those for plates.
Similarly, the normal strain components ε k can be expressed as follows and D k z = ∂ ∂z I 3 . Assuming a linearly elastic orthotropic material, the constitutive equations of the kth layer in the laminate reference coordinate system are written as is the vector of normal stresses, and the matrices of stiffness coefficients given bỹ are derived from those expressed in the layer reference system through a proper coordinate transformation [29]. Within the modeling framework proposed by Carrera [21], an entire class of 2-D shell theories can be employed by expressing the displacement vector u k through the Einstein notation as follows where ζ k = 2z k /h k is the dimensionless thickness coordinate of the layer (−1 ≤ ζ k ≤ 1), τ is the theory-related index, Fτ(ζ k ) are appropriate thickness functions defined locally for each layer, and is the vector of generalized kinematic coordinates in the assumed displacement model corresponding to index τ. Various ESL and LW theories of different order, which is a free parameter of the formulation, can be obtained by choosing the type of thickness functions and the range values of τ in Eq. (9). Both a class of LW theories and a class of ESL theories are implemented in this work. They are briefly presented in the following. A family of LW theories of variable order N is considered by assuming τ = t, b, r (r = 2, . . . , N) and selecting where Pr(ζ k ) is the Legendre polynomial of rth order. In so doing, the displacement variables u k b and u k t are the actual values at the bottom and top surfaces of layer k, respectively, and the interlaminar displacement continuity can be easily imposed as u k t = u k+1 b . Following the nomenclature suggested by Carrera [21], each member of this family is here denoted by the acronym LDN, which stands for (L)ayerwise (D)isplacement-based theory of order N. Note that the number of degrees of freedom for a LDN theory depends on the number of layers of the laminated panel and is given by The class of ESL theories considered here involves global thickness function as increasing power of the global thickness coordinate z, which is measured from the middle surface of the curved panel. Since the kinematics is layerindependent, the k index in Eq. (9) is dropped. Accordingly, t = 0, b = 1, and The related N-order member of the ESL family is indicated by EDN, which stands for (E)quivalent single-layer (D)isplacement-based theory of order N. Using Eq. (9), the geometric relations and constitutive equations of each layer of the curved panel are written, respectively, as and where the index s for the thickness functions in Eq. (14) has the same role of the index τ.
The equations of motion are here derived from the principle of virtual displacements (PVD), which, in absence of any external body and surface force, can be expressed as follows where Z k is the layer domain in the thickness direction. Using Eq. (9) and the geometric relations of Eqs. (13) into the above PVD statement, the dynamic equilibrium is expressed in weak form by the following equation is a thickness integral over the kth layer involving the product of thickness functions Fτ and Fs.
After integrating by parts the second term in Eq. (16), using the constitutive relations (14), and exploiting the arbitrariness of the virtual variation of each kinematic variable δu k τ over Ω k , the following set of differential equations in compact indicial notation is obtained: where L kτs and M kτs are 3 × 3 matrices of differential operators expressed as and It is noted that Eq. (18) is valid for any index τ and any layer k, and the summation for repeated index s is implied. Therefore, according to the order of the assumed shell theory and the number of layers of the curved panel, the expanded set of governing equations can be derived from Eq. (18) for each specific case under investigation. However, this set is not written in explicit form, since, according to Carrera's formulation, L kτs and M kτs are profitably viewed as invariant entities, which are directly expanded and assembled to build in a hierarchical and automatic manner the discretized stiffness and mass matrices of the free vibration problem, as explained later. For this reason, L kτs and M kτs are called the fundamental nuclei of the formulation [21].
Similarly, the arbitrariness of the virtual variation δu k τ over the boundary of the panel yields the following set of boundary conditions where B kτs is another invariant 3×3 matrix of differential operators called fundamental nucleus related to the boundary conditions. It is given for Neumann-type boundary conditions by where N k p and N k n are matrices containing the components nα and n β of the unit normal vector to the boundary as follows Since we are interested in the computation of the natural frequencies of the panel, a harmonic solution u k s =û k s e jωt is assumed. The corresponding differential eigenvalue problem is expressed by the equations and

Spectral collocation solution
A spectral collocation solution of the eigenvalue problem described by Eqs. (24) and (25) is sought. First, a grid of Chebyshev-Gauss-Lobatto (CGL) points (α i , β j ), i, j = 0, . . . , N CGL , over Ω k is introduced [23]. The discrete solution of the problem is assumed in the form of tensor product of one-dimensional expansions as follows whereû k αsmn ,û k βsmn andû k zsmn are the unknown values of the kinematic quantities at the grid points, i.e.,û k αsmn = u k αs (αm, βn),û k βsmn =û k βs (αm, βn) andû k zsmn =û k zs (αm, βn), and ψ l (l = 0, . . . , N CGL ) denotes the Lagrange interpolating polynomial relative to the given set of CGL nodes, so that ψm(α i ) = δ mi and ψn(β j ) = δ nj for i, j = 0, . . . , N. The approximation (û k αs ,û k βs ,û k zs ) is found by collocation; that is, (û k αs ,û k βs ,û k zs ) is required to satisfy the differential problem Eq. (24) along with the boundary conditions Eq. (25) at the CGL nodes. The derivatives in the differential operators L kτs and B kτs are approximated through the interpolation derivatives, corresponding to the evaluation at the grid points of the first-and second-order derivatives of the interpolants in Eq. (26) as shown below.
In order to simplify the writing and the coding process, the collocation equations are written in a matrix form. To this aim, the vectors U k αs , U k βs , U k zs ∈ R (NCGL+1) 2 of valuesû k αsmn ,û k βsmn andû k zsmn of the kinematic quantities at the CGL nodes are introduced as follows They are sorted so that increasing values of the coordinate α (β) of the curved panel correspond to increasing values of the index i (j). The values of the derivatives at the CGL points are also collected into vectors and are expressed in terms of the quantities in Eq. (27) by using the differentiation matrices D 1 and D 2 which contain the first-and second-order derivatives, respectively, of the interpolants evaluated at the CGL nodes [30]. We have D 2 = D 2 1 and the entries of D 1 are reported in [30]. For example, the evaluation of ∂û k αs /∂α at the collocation points is expressed as (D1 ⊗ I) U k αs , where I is the (N CGL + 1) × (N CGL + 1) identity matrix and ⊗ is the Kronecker product. Similarly, ∂ 2ûk zs /∂β 2 evaluated at the CGL nodes is given by (I ⊗ D 2) U k zs , and so on. By using the above notation, the following collocated equations expressed in the same indicial and invariant form of the original differential equations are obtained or, more compactly, are the (N CGL + 1) 2 × (N CGL + 1) 2 element matrices of the discrete stiffness-related fundamental nucleus K kτs arising from the spectral collocation process, and M kτs 11 = ρ k J kτs αβ (I ⊗ I) where e b is a (N CGL + 1) 2 boundary unit vector, i.e., a vector having all elements equal to zero except for one element equal to one corresponding to the panel edge under consideration. As outlined before, the final set of discretized equations involving the global stiffness and mass matrices of the panel and the final set of related boundary conditions are obtained directly from the discrete fundamental nuclei K kτs , M kτs and B kτs by a simple expansion and assemblylike procedure, which is graphically represented in Figure 2 for the stiffness nucleus K kτs . The same procedure is applied to the mass M kτs and the boundary B kτs nuclei. First, the discrete nuclei are expanded by varying the theory-related indices τ and s to yield and and M k and B k matrices are built in a similar way. Afterwards, the multilayered matrices are assembled by varying the index k. This yields and where K, M and B are simply summed layer-by-layer in the case of ESL theories, or they are assembled as shown in Figure 2 by enforcing the interlaminar continuity condition in the case of LW theories. Since the equations of motion must be enforced at the interior CGL points only and the boundary points must be explicitly defined in order to connect the differential problem with the set of boundary conditions, the final step of the spectral collocation process implies that Eqs. (36) and (37) are rewritten by partitioning the vector of unknown degrees of freedom U and the related matrices into quantities related to domain interior points (I) and quantities related to boundary points (B). Therefore, the previous eigenvalue problem can be expressed as From the second equation, we have Substituting Eq. (40) into Eq. (38), the natural frequencies of the curved panel can be obtained by solving the following generalized eigenvalue problem

Numerical examples
This section presents some illustrative numerical examples with the aim of showing the convergence properties, versatility and accuracy of the modeling and solution framework proposed in this work. The curved panels have side lengths a and b in the α and β direction, respectively, as shown in Figure 1.
In all examples, the boundary conditions of the panel are denoted by a four-letter notation numbered in a counterclockwise direction starting from edge α = 0. Both single-layer and multilayered panels made of isotropic and orthotropic material are considered. The properties for orthotropic material are the following: The mass density is taken as a reference dimensionless value of ρ = 1. For laminated panels, the layers are assumed to have the same thickness h k . The convergence rate of the spectral collocation solution is first examined in Figure 3, which reports the log plot of the absolute value of ∆ = λ NCGL − λ NCGL−2 , where λ NCGL denotes the fundamental non dimensional frequency parameter ω(a 2 /h) √︀ ρ/E 2 computed using N CGL colloca- tion points along α and β direction. N CGL is varied in the analysis from 6 to 20 with a step of 2. The computations are referred to a moderately thick (a/h = 10) square fully simply-supported (SSSS) spherical panel having [0/90/0] layup and two different shallowness ratios, a/Rα = 0.2 (moderately deep panel) and a/Rα = 0.02 (shallow panel). Values are obtained in both cases using the ED2 shell theory. It is shown that the spectral solutions exhibit a very fast convergence rate with an exponential decrease of ∆ as the number of CGL points increases. The difference is so small for N CGL ≥ 14 that it is swamped by the round-off errors of the numerical computations. The same spherical panel of the previous case but with fully clamped (CCCC) edge conditions is considered in Figure 4. The fundamental frequency of the problem is again computed using ED2 theory. It is observed that the clamped boundaries largely affect the convergence properties of the solution. Although converged results to five significant digits are obtained using a relatively low number of collocation points, clamped edge conditions slow down the rate of convergence. Apparently, no further improvement of the approximate solution occurs when N CGL > 14.
The adversing effect of clamped edges on the convergence properties of the present spectral solution is also found in other cases. For example, Figure 5 shows the difference between two successive approximations evaluated increasing by two the number of collocation points for the first four dimensionless frequencies λ i = ω i (a 2 /h) √︀ ρ/E 2 of a moderately thick (a/h = 10) square fully clamped [45/ − 45] 3 deep cylindrical panel with a/R β = 2. In this case, N CGL is varied from 6 to 18, and the results are obtained with a LD1 theory. It is noted that the convergence rate is similar to the case of a fully clamped spherical panel and is not significantly affected by the lamination lay-up, the shallowness ratio and the shell theory. Furthermore, no notable difference is observed as the number of mode to be estimated increases.
Similar trends of those shown in the previous figures have been observed in many other numerical cases, which are not presented here for the sake of brevity. It can be argued that sufficiently well converged values for the natural frequencies of the curved panels under study are obtained with N CGL = 14, which is the number of collocation points used in the following comparison analysis.  The accuracy of the modeling approach proposed in this paper is now evaluated and discussed by comparing the present solutions computed on models relying on shell theories of different topology and order with some reference solutions existing in the literature.
Two cases involving isotropic panels are first considered. One example deals with a SSSS square (a = b) cylindrical panel having two different span-to-thickness ratios a/h = 20 and a/h = 10 and three shallowness ratios a/R β = 0.5, 1 and 2. Table 1 lists the first five natural frequency parameters λ i = ω i (a 2 /h) √︀ ρ/E found using the present approach with ED2 and ED3 theories. Results are compared against those obtained by [15] using 3-D elasticity from a high-fidelity finite element model. It can be seen that the third-order shear and normal deformation theory predicts all five frequencies very close to 3-D results. Values from ED2 models are quite accurate for the fundamental mode. However, the accuracy deteriorates as the vibration mode number increases. The second case of isotropic panels considered here involves CCCC square spherical shells with a/Rα = 0.5 and length-to-thickness ratios a/h = 10 and a/h = 5. Non-dimensional frequency parameters λ i = ω i a √︁ ρ E computed using a fourth-order ED4 theory are shown in Table 2 and compared with the 3-D elasticity solutions obtained by [14] from Ritz models. A good agreement is observed even when a thick panel is considered.
Fully clamped composite laminated cylindrical panels are examined in Tables 3 and 4. Both cases are referred to a square moderately thick deep panel with a/h = 10 and a/R β = 2. Table 3 shows the first five natural frequency parameters for a cross-ply [0/90] 3 cylindrical panel as computed with ESL and LW theories of increasing order. The same analysis is carried out in Table 4 but for an angleply [45/ − 45] 3 panel. Present 2-D results are compared with those arising from a high-fidelity 3-D finite element analysis [19], which is assumed to be highly reliable. It is worth noting that ESL theories are rather inaccurate in both cases and for all frequencies under study, including the fundamental one. Increasing the order N of the ESL theory slightly improves the estimation, but the discrepancy between EDN and 3-D results is still significant even when N = 6 is adopted. The smallest difference in the results of EDN theories occurs in the fundamental frequency   parameter λ 1 for the cross-ply lay-up wherein ED6 predicts a frequency 5.2% higher than 3-D analysis. A dramatic increase in accuracy is observed when 2-D layerwise theories are implemented. This is due to the capability by LW models of correctly predicting the through-the-thickness zig-zag behavior of the displacements in correspondence of each layer interface. The accuracy of LD1 is already acceptable with differences of about 1% with respect to 3-D models. If an improved accuracy is required, more costly higher-order LDN theories can be used. Table 5 shows the same results as in the Table 3 for a CSCS cylindrical panel. A subset of ESL and LW models that can be implemented in the present modeling framework is selected. It can be seen again that the accuracy of EDN theories cannot be significantly improved beyond a certain limit by increasing the order N. Instead, LDN theories are capable of providing frequency parameters very close to 3-D values. As before, the degree of accuracy improves by using LW models of higher order.
In the last numerical example, a simply-supported SSSS square spherical panel is considered. It is composed of three layers oriented at [0/90/0]. The analysis is carried out for a/h = 10 and a/h = 100 (thin panel) and for different values of a/Rα including the plate case when Rα → ∞. The fundamental frequency parameter λ = ω(a 2 /h) √︁ ρ E2  is reported in Table 6 for each case. Comparison is provided with values from [24] where a layerwise approach is adopted in combination with a multiquadric radial basis function method. Results are in good agreement for both thin and thick case. In particular, it is shown that ESL models give very accurate predictions in the thin case without the need of relying on costly LW models.

Conclusions
An advanced modeling and solution framework for the prediction of natural frequencies of both thin and thick,  deep and shallow cylindrical and spherical panels with different boundary conditions and arbitrary lamination lay-up is presented. The approach combines the capability of building hierarchically 2-D shell models based on equivalent single-layer and layerwise theories of variable order and the efficiency and simplicity of the spectral collocation solution method. The formulation is expressed in a compact matrix form which can be directly and easily coded in any modern mathematical software package. It has been shown through comparison with a fully 3-D analysis that stringent requirements on the accuracy of the computed frequency values can be satisfied only by 2-D high-order layerwise models, in particular when thick and deep anisotropic curved panels are considered. However, the present approach can also allow to build less costly low-order equivalent single-layer models when a more refined analysis is not required for the specific case under investigation.