Application of the ps−Version of the Finite Element Method to the Analysis of Laminated Shells

The development of accurate and efficient numerical methods is of crucial importance for the analysis and design of composite structures. This is even more true in the presence of variable stiffness (VS) configurations, where intricate load paths can be responsible for complex and localized stress profiles. In this work, we present the ps−version of the finite elements method (ps−FEM), a novel FE approach which can perform global/local analysis through different refinement strategies efficiently and easily. Within this framework, the global behavior is captured through a p−refinement by increasing the polynomial order of the elements. For the local one, a mesh−superposition technique, called s−refinement, is used to improve locally the solution by defining a local/fine mesh overlaid to the global/coarse one. The combination of p− and s−refinements enables us to achieve excellent accuracy−to−cost ratios. This paper aims to present the numerical formulation and the implementation aspects of this novel approach to VS composite shell analysis. Numerical tests are reported to illustrate the potential of the method. The results provide a clear insight of its potential to guarantee fast convergence and easy mesh refinement where needed.


Introduction
The anisotropic nature of composite laminates is responsible for structural responses characterized by complex couplings between extension, bending, torsion, and shearing modes [1]. Such effects are particularly relevant when laminates with curvilinear reinforcements [2][3][4][5] are of concern.
These effects require appropriate numerical tools to be available to conduct accurate and reliable analyses. This is even more true as composite thin shells are increasingly used in primary load−carrying structural components. Hence, the interest of the designers is not restricted to the global response, but includes also local effects, such as stress profiles in critical regions due to geometric discontinuities, material interphases, and areas affected by load introduction. In this framework, the finite element method (FEM) is the most commonly used tool for the analysis and design of composite structures. The main advantages can be found in the high flexibility of modeling complex domains, the possibility of considering arbitrary material properties, and handling different loading/boundary conditions. Over the years, different strategies have been developed to maximize the accuracy of the FE solution, while keeping the computational cost at a minimum. The classical approaches to improve accuracy are represented by the h− and p−refinements [6]. In the former, accuracy is increased by reducing the mesh size h, while the interpolation order of the elements p is fixed [7]. In the latter, the polynomial degree p is increased, while the mesh resolution h is left unchanged [8]. The h− and p−refinements provide two systematic ways to improve the precision of the FE solution and yield excellent results when global quantities are of interest. Different applications of these strategies to advanced composite shells and plates can be found in the literature, such as sandwich plates [9], functionally The outline of the paper is as follows: the shell mathematical model is presented in Section 2; Section 3 discusses the formulation and implementation aspects of the proposed FEM; in Section 4, the method is applied to several test cases involving laminated shells. Finally, the main findings of this study are summarized in Section 5.

Theoretical Model
This section illustrates the theoretical framework for the analysis of laminated composite shells. In the first part, the shell geometry is discussed by presenting the relevant equations employed next in the finite element approximation. Then, the constitutive law for handling the case of curvilinear fiber paths is presented. Finally, the shell governing equations are derived in the context of a total Lagrangian description.

Geometry Description
The geometric description of shells is developed under the theory of surfaces, which provides a general framework to define panels with arbitrary configurations and curvatures.
The reference surface Ω is assumed to be coinciding with the shell's midsurface and is parametrized by two arc−length coordinates ξ 1 ∈ [0, a] and ξ 2 ∈ [0, b], where a and b are the curvilinear lengths of the shell edges. The normal direction to Ω is parametrized by the coordinate ζ ∈ [−t/2, +t/2], where t is the thickness of the shell. It is assumed that the triad (ξ 1 , ξ 2 , ζ) forms an orthogonal curvilinear coordinate system such that ξ 1 and ξ 2 are aligned with the principal lines of curvature on the midsurface, as shown in Figure 1.
Considering a global reference system (O, x, y, z), any point on Ω can be described by the position vector r: r(ξ 1 , ξ 2 ) = r x (ξ 1 , ξ 2 )e x + r y (ξ 1 , ξ 2 )e y + r z (ξ 1 , ξ 2 )e z (1) where r x , r y and r z are definite, continuous, single-valued functions of ξ 1 and ξ 2 representing the components of r on the global reference system with unit directors e x , e y and e z . The tangent vectors on Ω to the two coordinate lines, ξ 1 and ξ 2 , are obtained by differentiation of Equation (1): where the comma denotes differentiation with respect to the coordinate ξ α . Accordingly, the unit normal to the reference surface is defined as: where × denotes the cross product and | · | the Euclidean norm. The unit normal in Equation (3) is assumed to be common to all the shell laminas Ω (ζ) , see Figure 1.
Considering a lamina at a distance ζ from Ω, the position vector of the generic point on Ω (ζ) is expressed as: r (ζ) (ξ 1 , ξ 2 , ζ) = r(ξ 1 , ξ 2 ) + ζn(ξ 1 , ξ 2 ) (4) while the corresponding tangent vectors along the curvilinear lines ξ 1 and ξ 2 are: ,α = r ,α + ζn ,α for α = 1, 2 By application of the Weingarten-Gauss relations [1]: the following relation between the tangent vectors on Ω (ζ) and Ω holds: where R α are the principal radii of curvature of the reference surface and are defined as [44]: where · denotes the dot product. In general, the two principal radii of curvatures R 1 and R 2 are not constant and can be a function of the coordinate lines ξ 1 and ξ 2 .
Referring to Equation (4), the differential position vector is given as: From Equation (9), the square distance dl between two arbitrary points in the shell volume P(ξ 1 , ξ 2 , ζ) and P (ξ 1 + dξ 1 , ξ 2 + dξ 2 , ζ + dζ) is computed as: where a (ζ) α are the Lamè parameters defined as: with: For the development of the shell mathematical model, the geometric measures are expressed in terms of reference surface Ω quantities. Accordingly, the element lines and cross section areas along the coordinates ξ 1 and ξ 2 are given as: while the element surface and volume are given by: The geometric formulation presented here furnishes the background to develop the finite elements discussed next, where plates and shells with variable curvature can be taken into account.

Curvilinear Fiber Path Description
The shell elements considered here can be applied to the analysis of isotropic and composite structures. For generality purposes, the model is developed to allow non−uniform elastic properties to be considered along the surface. Hence, new configurations such as variable stiffness (VS) laminates can be studied within the proposed framework. Specifically, a laminate with VS properties is achieved by stacking together plies with curvilinear reinforcing fibers [2]. Different strategies have been proposed in the literature for describing the path of the fibers. Without loss of generality, the approach considered here relies on the Lagrange polynomials introduced in [4]. According to this mathematical formulation, the curvilinear paths are described by a set of parameters T mn according to: (15), it is possible to understand the parameters T mn as the fiber orientation angle at the control points, as one can see that θ ξ = T mn . The information of fiber orientation on a the generic VS ply can thus be defined through a matrix T ∈ R M×N . This enables us to specify VS lamination sequences very concisely as [T 1 /T 2 / . . . /T N p ], where N p stands for the number of plies.
An alternative formulation is also possible where the control points are defined at a quarter of the laminate and the resulting fiber orientation distribution is reflected to the rest of the domain. In this case, the expression of Equation (15) modifies as: where the new coordinates areξ 1 = ξ 1 − a/2 andξ 2 = ξ 2 − b/2, while the control points

Shell Description
The displacements of an arbitrary point on the shell domain Ω are denoted by u 1 , u 2 and u 3 . These components correspond to the displacements along the ξ 1 , ξ 2 and ζ directions, respectively. The component u 3 is taken positive in the outward direction from the center of the smallest radius of curvature. Considering only von Kàrmàn nonlinearities, the components of the Green-Lagrange strain tensor in the orthogonal curvilinear coordinate system (ξ 1 , ξ 2 , ζ) are [45]: where αα and γ αβ are the normal and shearing strains, respectively. The shell kinematics is described according to the first shear deformation theory (FSDT), so: where u, v and w are the displacement components of an arbitrary point on Ω along ξ 1 , ξ 2 and ζ, respectively, while φ 1 and φ 2 are the rotations. The choice of FSDT is due to the easier FE implementation compared with the classical Kirchhoff theory (CLT). At the same time, FSDT allows shear deformability to be taken into account with a small number of unknown fields. Higher−order theories are not considered here, but can be of practical interest to guarantee improved description of the thickness−wise response, while avoiding the use of the shear factor. The strains of the kinematic model at hand can be obtained upon the substitution of Equation ( 13 γ 0 23 } T is available in the Appendix A. The membrane resultant on the cross section dA 2 , i.e., the section normal to ξ 1 , is: where the definition of the resultant N 11 follows from the equation above. The other resultants are derived in a similar manner, leading to: where K s is the shear correction factor, which is taken as equal to 5/6. The shear stress resultants, N 12 and N 21 , and the twisting moments, M 12 and M 21 , are different due to the curvature terms. For shallow and moderately thick shells, i.e., h/R α < 1/20, the simplifying assumption N 12 ≈ N 21 and M 12 ≈ M 21 can be introduced. Based on the assumption above, the constitutive law reads: The equations of motion are derived by referring to Hamilton's principle: where K, U and V are the kinetic and elastic energies and the potential of the applied loads, respectively. For the shell model developed in this work, these energy quantities are defined as follows: where the dot denotes the time derivative, q + is the pressure applied on the upper surface (ζ = h/2), ∂Ω 1 and ∂Ω 2 are the boundaries at ξ 1 = const and ξ 2 = const, respectively. The terms I α are the moment of inertia defined as: where ρ is the mass density.

The ps−Version of the Finite Element Method
The finite element method (FEM) developed in this work is a combination of two FE schemes, i.e., the p−version of the finite element method (p−FEM) and the s−version (s−FEM). The former is an extension of the conventional FEM, or h−version (h−FEM), where the accuracy is increased by increasing the interpolation order (p−refinement). The latter is an FE scheme where arbitrary local improvements of the mesh resolution are possible through the adoption of advanced mesh superposition techniques (s−refinement). The two features of these numerical approaches are put together into an extended FE framework, called the ps−version of the finite element method (ps−FEM), where simultaneous p− and s−refinement (ps−refinement) can be performed to adaptively adjust the interpolation order p and the size h of the elements. Within the ps−FEM, the construction of the polynomial space S p is inspired from the p−FEM, while the design of the mesh ∆ follows the ideas implemented in the s−FEM, see Figure 2. In the following, the two main pillars of the ps−FEM, i.e., p−FEM and s−FEM, are presented. Then, the ps−FEM is introduced and applied to composite shell problems by referring to the equilibrium equations of Section 2.3.

p−Refinement
The core idea of the p−refinement strategy is to progressively increase the polynomial order p of the shape functions to improve the quality of the numerical solution. One crucial aspect regards the worsening of the conditionality of the global stiffness matrix as p is increased. In the p−FEM, this issue can be mitigated by constructing the set of shape functions starting from Legendre polynomials. The properties of these functions enable better numerical conditioning in comparison to the Lagrange polynomials commonly used in the h−FEM. High−order polynomials can be considered in the p−FEM-even beyond p = 3-, without encountering numerical issues.
Considering a generic two−dimensional solution field φ depending on x and y, the FE discretization in the elemental domain Ω (k) is: where −1 < {ξ, η} < 1 are the nondimensional coordinates in the computational domain Ω st , χ (k) is an array of functions mapping the standard element in the k−th element of the mesh ∆, see Figure 3, the coefficients p and c ij are the element interpolation order and degrees of freedom, respectively, while f i and f j are polynomial shape functions. In the p−FEM the polynomial space S p is constructed as: where P n (ξ) is the Legendre polynomial functions of order n: From Equation (27), it is possible to distinguish two families of functions describing different deformation modes, i.e., nodal and internal ones. Nodal modes are represented by f 1 (ξ) and f 2 (ξ), which are Lagrange linear interpolation polynomials. Internal modes are reproduced by the higher−order terms, f α (ξ) (for α ≥ 3), which are based on the integrals of Legendre polynomials. The combination of nodal and internal modes leads to different types of two−dimensional shape functions, i.e., nodal (two nodal modes), edge (one nodal mode and one internal mode), and face (two internal modes) shape functions, see Figure 4.
The series expansion in Equation (27) offers several advantages over the classical Lagrangian polynomial base. The first one is that it inherits the orthogonal properties of the original Legendre polynomials, indeed:  This property provides the stiffness matrix with a quasi−diagonal structure, which significantly eases the solution process and enables the use of shape functions with a very high polynomial order. The second feature regards the hierarchical nature of these functions: p−refinements implies adding new higher−order terms in Equation (27), while the low−order ones remain unchanged. As a consequence, the final stiffness matrix is obtained by successively adding rows and columns associated with new levels of p−refinements.
Relatively few elements are in most cases sufficient to meet the desired levels of accuracy. At the same time, the method should be capable of appropriately representing the geometry of the structure with a limited number of elements. For this purpose, blending functions [46,47] are usually combined with the p−FEM. The mapping procedure for the transformation of the standard element in Ω st to the k−th element of the mesh in Ω (k) is: where Q (k) x (ξ, η) and Q (k) y (ξ, η) are general nonlinear mapping functions, (X α , Y α ) are the nodal coordinates of the element, while e αx and e αy with α = 1, . . . , 4 are the functions defining the curves of the element edge, as illustrated in Figure 3.
From the mapping of Equation (30), the derivatives with respect to the physical coordinates x and y are computed as: αβ are the components of the Jacobian matrix J (k) . Line integrals at the four sides of the element read: while the surface ones are: where F is the generic integrand, F is obtained from F by replacing x and y with the mapping functions in Equation (30), while J = detJ (k) is the determinant of the Jacobian matrix.

s−Refinement
The main problem in performing local mesh refinements relates to the need to generate a transition between refined and unrefined regions. Different approaches are available in the literature to address this issue, for example transition elements or multi−point constraints approaches [24][25][26][27]. The s−refinement strategy offers the advantage of simplifying this process by allowing an element size reduction in the desired regions only. This result is achieved through the definition of an independent local/fine mesh which is superimposed to a global/coarse one [30].
Using this idea, the final FE approximation φ can be represented as: where φ G is the global mesh solution defined in Ω, while φ L is the local mesh solution defined in Ω L ⊂ Ω, see Figure 5. Note that the mesh superposition technique allows for incompatible discretization between global and local mesh. This gives an extremely high level of flexibility when performing local h−refinements, as no transition regions [24][25][26] or multi-point constraints [27] are required. This concept of solution superposition can be extended to multi−level refinements [41]. In this case, the elements of the local mesh can be further refined by superposing one over the other multiple levels of overlaid meshes. In this case, the final FE solution φ becomes: where φ Two conditions, i.e., compatibility of the basis functions and their linear independency, are required to apply this multi−level decomposition of the solution field φ.
The first condition implies C 1 −continuity within each elements and C 0 −continuity across the element boundaries. The C 1 −continuity is satisfied by construction. On the contrary, the inter−element continuity is not guaranteed and needs to the imposed. This is achieved by enforcing homogeneous Dirichlet boundary conditions on the boundary of the overlaid meshes, as depicted in Figure 6a. The second condition on the linear independency is required to avoid singularities in the stiffness matrix. In general, the redundant degrees−of−freedom can be removed during the factorization process by elimination of the equations with zero pivots [30]. If the p−FEM is employed for the global and local meshes, this is avoided by ensuring that shape functions of the same type (nodal, side, face) and polynomial order p appear only once in regions with multiple meshes. This idea is graphically illustrated in Figure 6b.

ps−Refinement
The ps−refinement procedure is a combination of the p− and s−FEM approaches. In this framework, the order of interpolation p can be increased with no numerical issues owing to the properties of Legendre polynomials; at the same time, the elements size h can be adaptively reduced by overlaying different levels of superposition meshes s.
In this work, the ps−refinement is used to develop an advanced FEM, called the ps−FEM, which is exploited to perform global/local analysis of laminated shells. In particular, the p−refinement strategy is exploited to capture the smooth features of the solution at global scale, such as the global deformation field, while the s−refinement approach is employed to capture local effects, such as stress concentrations.
The refinement process is outlined in Figure 7. First, a global/base mesh ∆ G is defined. The resolution of ∆ G is chosen to guarantee appropriate definition of the loading and boundary conditions, and to avoid the description of the geometry using distorted elements. Then, a p−refinement is performed until the required accuracy on the global response is reached. Finally, the solution is s−refined with local meshes in the regions where enhanced interpolation capability is required. The ps−FEM shares similar convergence features with the hp−FEM. Indeed, by increasing the order of the polynomial expansion in combination with overlaid meshes, it is possible achieve an exponential decaying global approximation error. This is true even when the solution presents steep gradients or singular points. More insights on convergence properties can be found in the work of [42].
In the proposed framework, three different combinations of p− and s−refinement are allowed, namely linear, uniform, and graded ps−refinements. In the linear case [39], the global solution is represented by a coarse high−order mesh, while locally a multi−level s−refinement is performed with low−order meshes. The uniform and graded ps−refinement strategies [41] are an extension of the previous one, where high−order local meshes are employed. In the first case, the local meshes have the same order p of the global one, while in the second case p is different from level to level.
In principle, s−refinement can be carried out considering local meshes with arbitrary resolution and orientation, as originally done in [30]. A simplified approach is implemented here, inspired by the work of [41]. In particular, the elements of the refined region are divided to half of the original size. These elements are then used to define the superposed mesh at the upper level. This procedure allows for a quicker refinement, with no need to iterate to correlate the generic coordinates of the local and global element domain [30].
In the ps−FEM, the shell generalized displacements for the k−th element with interpolation order p at the s−level mesh are: where p α and c (α) (α = 1, 2, 3, 4, 5) are the element expansion order and unknown amplitude coefficients for the displacement fields, while χ (k,s) is the vector collecting the mapping functions between the element computational space Ω (k,s) st and element physical space Ω (k,s) . This mapping procedure is illustrated in Figure 8. For the global mesh elements (s=0), the vector of mapping functions is defined by Equation (30). For the superimposed mesh elements (s > 0), the following sequence of mapping is operated: where Ψ (k,s) is the vector collecting the mapping functions between the computational domain of the element and the underlying one, as shown in Figure 8. Each element is therefore provided with two set of mapping functions, the global one χ (k,s) and the local one Ψ (k,s) . These functions are required for defining the derivatives and integrating at element level. It is noted that the ps−FEM can be applied for the numerical solution of any mathematical model, such as higher−order shear deformation theories, as well as three−dimensional elasticity theories. The effectiveness of FE techniques relying on mesh superposition has been discussed in the literature for such problems. Examples can be found in [42], where the hp−d−FEM is employed to solve 3D linear elastodynamic problems, and in [32,35] where the s−FEM is employed in the context of multiple model methods.
Using the approximation of Equation (36), the energy contributions of Equation (24) are: where use is made of the index notation presented in [48]. According to this notation, a pq represents a vector, A (pq)(rs) , A (pq)(rs)(mn) , A (pq)(rs)(mn)(ij) are second−, third− and fourth−order fourth−order arrays. By application of Hamilton's principle to the generic element, the following kernel is obtained: where u (k,s) is the vector collecting the element degrees of freedom, M (k,s) , K (k,s) , P (k,s) are the element stiffness, mass matrices, and load vector. Upon assembly of the contributions of all elements, the governing equations are obtained as follows: where M, K and P are the global mass, linear stiffness matrices and vector of external loads, respectively, while N 2 and N 3 are the nonlinear contributions due to cubic and quartic terms in the elastic energy, respectively. The set of ODEs reported in Equation (40) describe the nonlinear dynamics of the shell and can be solved via implicit or explicit integration techniques. In this work, the discussion is restricted to the quasi−static case, so the relevant equations simplify to: for which the solution is computed by referring to the Newton-Raphson method.

Results
This section presents applications of ps−FEM to laminated shell problems. The method has been implemented in a Matlab environment, while all analysis have been run on a laptop with the following characteristics: 1.4 GHz Intel Core i7 processor, 16 GB 1867 MHz LPDDR3 memory, Intel HD Graphics 6,151,536 MB.
In the first part of the section, a series of validation studies is proposed. Then, the numerical technique proposed herein is exploited to solve exemplary test cases for laminated shell problems. The potential of the ps−refinement is illustrated against classical h− and p−refinements for linear and geometrically nonlinear problems.

Test Case 1: Vibrations of Elliptical Shells
The first study aims at validating the shell model implemented and, for this purpose, the vibration response is studied for three elliptical cylindrical shells. The test case is taken from the literature [49] and covers the case of shells characterized by different eccentricities e. The position vector of the reference surface is expressed as: r(ξ 1 , ξ 2 ) = ξ 1 e x + R A sin (ξ 2 )e y + R B cos (ξ 2 )e z for ξ 1 ∈ [0, L] and ξ 2 ∈ [0, 2π] (42) where R A and R B are the semi−major and semi−minor axis, respectively, while L is the length. A sketch of the shell is provided in Figure 9, while the geometric configurations considered in this study are summarized in Table 1. The shells are made of an aluminum whose properties are E = 68, 950 MPa, ν = 0.3, and ρ = 2766 kgm/ 3 . Clamped−free boundary conditions are considered.  The numerical model consists of a coarse mesh with five elements along the circumferential direction and two along the axial one, see Figure 10. The polynomial order is p = 10. The results are reported in Table 2 for the first ten natural frequencies. A comparison is illustrated against the reference results taken from [49], where the same shell kinematics and geometric formulation were adopted. In [49] the shell governing equations are expressed in strong−form and are solved using the generalized differential quadrature method. Excellent matching is achieved for the reported frequencies, where the error is evaluated as E % = |ω − ω ref |/ω ref × 100. For the reported modes, the percentage errors are of the order of E % = 10 −2 − 10 −3 . The almost perfect agreement gives evidence of the correct implementation of the shell geometry and kinematics.

Test Case 2: Vibration of Variable Stiffness Plate
The second study aims at validating the ability of the proposed numerical code to handle the case of structures with non−uniform stiffness. In this regards, the example presented herein provides a validation of the correctness of the constitutive law implemented [50].
A free vibration analysis is conducted by considering a variable stiffness (VS) rectangular laminate, as illustrated in Figure 11. In particular, the geometry is defined by a = b = 300 mm and t = 1.2 mm. The material is characterized by the following elastic properties: E 11 = 10,000 MPa, E 22 = 9000 MPa, G 12 = G 13 = G 23 = 5000 MPa, ν 12 = ν 13 = ν 23 where different values of T 12 define different fiber paths. Starting from the configurations defined by Equation (43), the local fiber orientation is obtained by referring to Equation (16). The VS plate is clamped at the four edges and is modeled with one single element of polynomial order p. The results are presented in terms of the first nondimensional frequency Ω = ωa ρ/E 22 in Table 3. The comparison against the reference results is provided by considering different degrees of p−refinement. As shown, the FE results converge quickly upon p−refinement with frequencies converging from above. This behavior is motivated by the hierarchical properties of the Legendre expansion [6]. It is important to note that different FE discretizations may imply different distributions of the elastic properties leading to nonmonothonic convergence of the solution. The comparison of the converged frequencies with the reference ones [50] provides differences below 1%, thus illustrating the possibility of successfully handling structures with non−uniform elastic properties within the present framework.

Test Case 3: Static Analysis of Plate with Cutout
This test case deals with the static analysis of plate with circular cutout loaded in traction. The presence of the cutout determines the onset of stress concentrations, thus providing a useful mean for investigating the validity of the proposed ps−refinement strategies.
The plate, an illustration of which is reported in Figure 12, is characterized by dimensions a = 75 mm, b = 50 mm and r = 10 mm. The elastic properties considered for the material are: E 11 = 1,500,00 MPa, E 22 = 9000 MPa, G 12 = G 13 = G 23 = 5000 MPa, ν 12 = ν 13 = ν 23 = 0.32. A quasi−isotropic layup [±45, 90, 0] s is considered and each ply has thickness t ply = 1 mm. The traction load is prescribed at two parallel edges with magnitude N xx = 100 N/mm. The FE model is restricted to a quarter of the full structure owing to the double symmetry of the problem. The s−refinement strategy is considered for designing the mesh of the numerical model. In particular, a base mesh having five elements with polynomial order p is superposed with s = 5 levels of overlaying meshes of decreasing size, as illustrated in Figure 13. In the validation, the order p is increased until convergence of the global solution is reached. The number of overlaying layers s is selected to have a proper representation of the local response, while the resolution of the base mesh (number/size of elements) is chosen to guarantee the appropriate subdivision of the domain for performing s−refinements in the regions of interest. Specifically, the overlaying meshes are placed around the corners of the cutout, where stress gradients are expected. Three FE models employing the ps−refinement strategies illustrated in Figure 7 are set up for the validation study. The interpolation order of the superimposing elements are p s = 1 and p s = p for the linear and uniform ps−refinement, respectively. Regarding the graded ps−refinement, the polynomial order p s is set according to the following law: where m is a constant that controls the rate of decrease of p for increasing levels s. The law of Equation (44) will be adopted in the rest of this work for the graded ps−refinement strategy, unless otherwise specified. In the present study, the ps−refinements are performed fixing the mesh and increasing the polynomial degree of the elements on the global mesh p, while the order of the overlaying elements p s is adjusted according to the specific refinement strategy, i.e., linear, uniform, or graded. The validation is carried out using a FE model realized in Abaqus with S4 shell elements. A preliminary convergence study was carried out with this model. Based on this study, the total number of elements has been taken equal to 350,000 elements (h−refinement) to guarantee convergence. The corresponding number of membrane degrees of freedom is 7 × 10 5 , approximately. Results are shown in Table 4 for different types and levels of ps−refinements, and are presented in terms of nondimensional elastic energy: From Table 4, it can be seen that the nondimensional elastic energies of the three FE models approach the reference one as p is increased. In particular, one can note a quick monotonic convergence from below, i.e., the strain energy increases as the FE models are more refined. This is an expected behaviour as the numerical model becomes more flexible with an increasing number of degrees of freedom. For the linear and graded ps−refinement strategies, convergence is reached with p = 5. Regarding the uniform ps−refinement, convergence is achieved with a slightly lower polynomial order (p = 4). This is due to the use of high−order elements in the upper mesh levels which enable us to obtain a higher approximation capability. In Figure 14, the convergence curves are presented for the stress resultant N xx measured at (ξ 1 , ξ 2 ) = (0, r) for the three FE models considered for the validation. In all cases, one can note a strong and non−monotonic convergence of N xx to the reference value of 378.26 MPa. The curves in Figure 14 provide insights on the local convergence rate of the ps−refinement strategies. For the present problem, convergence of the stress measures is reached for approximately 10 3 degrees of freedom, which is much lower than the one required by the reference model where a h−refinement was employed. The consistency of these results, both in terms of global and local quantities, provide a validation for the ps−refinement strategies implemented in this work.

Example 1: Vibration and Buckling of a Highly Anisotropic Laminated Plate
The first application of the ps−FEM is an example taken from [51], and relates to the vibration and buckling analysis of a highly anisotropic composite plate. The following data are considered for the analysis: a/b = 1 and a/t = 100,000, see Figure 11. The material has elastic properties given by: E 11 = 393,000 MPa, E 22 = 5030 MPa, G 12 = G 13 = G 23 = 5240 MPa, ν 12 = ν 13 = ν 23 = 0.31, ρ = 1500 kg/m 3 . The plate is a laminate with a single ply oriented at θ = 45. Simply supported boundary conditions are considered along the four edges, while uniform edge−shortening conditions are considered for the buckling analysis.
This problem has been studied in previous efforts in the literature due to its challenging convergence features. The eigenvalues of the problem, both for free vibrations and buckling, tend to convergence with a slow rate as a consequence of the high orthotropy ratio E 11 /E 22 = 73.36, whose effects are exacerbated by the 45−oriented ply. Indeed, this orientation promotes strong bending/twisting elastic couplings with highly localized gradients in the modal shapes.
Two meshes are realized for conducting the numerical tests. The first one is a mesh with 5 × 5 quadrilateral elements, see Figure 15a. The second one is a s−refined mesh, whose features are outlined in Figure 15b; the s−refinement is obtained starting from the base mesh 5 × 5 and adding 15 layers of overlaying meshes. The refinement is conducted to allow local effects arising from vibration and buckling modes to be accurately captured. Starting from these two meshes, five FE models are developed to illustrate the effect of different refinement strategies. The first two models, h− and p−models, employ the mesh in Figure 15a and implement the h− and p−refinement strategies, respectively. The other three models, ps−L−, ps−U−, and ps−G−models, are constructed from the s−refined mesh in Figure 15b and adopt the linear (L), uniform (U), and graded (G) ps−refinement strategies, respectively. A summary of the FE models is reported in Table 5, whose nomenclature will be used for the other examples reported next. In absence of available exact solutions, the first step consists of obtaining highly accurate results to be used as a reference for comparing the five models later. The ps−U−model is used for this scope. The results are summarized in Table 6 in terms of nondimensional frequency and buckling load, defined as: (46) where N xx is the force per unit length on the loaded edges. As shown, a value of p = 7 guarantees convergence up to the third digit for both vibration and buckling parameters. These values are retained as references for computing the errors obtained with different refinement strategies.
The vibration and buckling mode shapes are illustrated in Figures 16 and 17. Both modes are characterized by one single−skew half−wave, where stretching occurs in the direction of fibers, i.e., θ = 45. From Figures 16a and 17a, one can note that the deflected shape w is smooth in most of the domain apart from the two corners. Here, it is possible to note highly localized effects, as evident from the plots of the twisting moments M xy in Figures 16b and 17b where strong stress concentrations are observed. This will have drastic consequences on the convergence of the fundamental frequency and critical buckling load.  The relative errors percentage from the numerical models of Table 5 are summarized in Figure 18, where an additional curve in the results of [51], i.e., Ritz p−refinement, is reported. The plots report the logarithmic relative errors percentage evaluated against the reference solutions. Despite quantitative distinctions, similar trends are observed for the frequency and the buckling results, the latter requiring more degrees of freedom for a fixed corresponding error. From Figure 18, the fast convergence rate achieved with the uniform and graded ps−refinements is clearly visible. In particular, the ps−U− and ps−G−refinements outperform FE models implementing standard h− and p−refinement approaches. This effect is clear from the error−to−degrees−of−freedom ratio of the nondimensional frequency and buckling load in Figure 18, with the ps−U− and ps−G−models presenting the curves with the steepest slopes. The excellent convergence properties just shown are made possible by the flexibility of the approach presented here: the degrees of freedom can be employed to refine specific areas of the domain with no need of generating transition regions, where refinement is not useful in terms of solution accuracy, but is a mean to generate compatible meshes.
The plots of Figure 18 allow us to distinguish regions where the adoption of a refinement strategy is advantageous over the others. For instance, percentage errors E % of the order of 0.1-1% suggest the adoption of a p−refinement strategy. On the contrary, if stricter requirements are set, the ps−refinements provide the best mean for solving the problem, with uniform and graded strategies as the most effective ones.

Benchmark results
The previous section demonstrated the advantages of the ps−FEM in deriving highly accurate solutions with relatively low computational costs. The proposed ps−FEM framework is then exploited to provide reference results to be used for comparison purposes in future studies in the field.
The analyses are performed using the ps−U−model considered in Table 6. The mesh is illustrated in Figure 15b, and the corresponding total number of degrees of freedom is approximately 2 × 10 4 .
A summary of the results is provided in Table 7 in terms of the nondimensional eigenvalues according to Equation (46). The results are compared with the predictions obtained using the Ritz method proposed in [51]. Specifically, the Ritz predictions are obtained by expanding the unknowns with 250 × 250 Ritz functions, corresponding to a total of 6 × 10 4 degrees of freedom, i.e., six times higher than the present model. As shown, the eigenvalues obtained via ps−U−model are always smaller than the Ritz ones. Owing to the convergence properties of the solution (the eigenvalues converge from above) the present results are thus more accurate. The maximum differences are reached for the plate with E 11 /E 22 = 73.36 and θ = 45, which corresponds to the largest degree of anisotropy. Note, for the buckling load, the difference is as high as 0.40%, approximately, with respect to an already refined solution. The difference between the present and reference results becomes smaller as the degree of anisotropy is reduced either by reducing the off−axis angle or the orthotropy ratio. The response of a cylindrical laminated shell is investigated under the effect of a static point load. The shell has an elliptical cross−section with semi−axes R A = 1000 mm and R B = 500 mm, total thickness t = 10 mm, and length L = 2000 mm. The sketch of the structure is reported in Figure 9. The material has elastic properties E 11 = 150,000 MPa, E 22 = 9000 MPa, G 12 = G 13 = G 23 = 5000 MPa, ν 12 = ν 13 = ν 23 = 0.32, while the stacking sequence is given as: [±T] s with T = [30,45,30] where the fiber path is obtained via Equation (15). The shell is clamped at both ends and is loaded with a transverse point load with magnitude P = 100 N and applied at (ξ 1 , ξ 2 ) = (L/2, 3/2π). This loading condition is of interest due to the localized phenomena associated with the concentrated force and finds practical application in tests such as the single perturbation load analysis (SPLA) [52], used for the assessment of shell imperfection sensitivity. Two FE models are developed and compared. The two models share the same number of degrees of freedom, 10 4 , approximately, but have different features for the meshing strategy. The first model relies upon a p−refinement strategy (p−model) and is presented in Figure 19a. Its mesh has 8 × 8 elements, i.e., four elements along the circumferential and the axial directions, with a polynomial order p = 11.
The second model adopts the s−refined mesh illustrated in Figure 19b. The base mesh is unaltered with respect to the p−model. A local refinement is introduced with s = 5 layers of overlaying meshes in the surroundings of the loaded region. The uniform ps−refinement strategy (ps−U−model) is considered, with the polynomial order of the base elements set to p = 5.   Figure 20b, where the discrepancies between the two models are noticeable. An excellent prediction of the stress gradient, which is particularly exacerbated close in the loaded region, is achieved via ps−refinement. The comparison between the peak values of the two models reveals an underprediction of the p−model that is close to 60%.
The contour of the stress is reported in Figure 21. As shown, the p−model solution presents severe oscillations in the surroundings of the stress peak (Figure 21a). This phenomenon is typically observed in high−order approximations when the solution displays steep gradients or discontinuities. Oscillations are greatly attenuated in the case of the ps−refined solution, as revealed by Figure 21b. These results provide a clear insight into the advantages of the proposed strategy. By recalling that the two models rely upon the same number of unknowns, evidence is given on the superiority of the ps−refinement strategy for problem characterized by strong gradients. As a matter of fact, the combination of p− and s−refinements allows us to combine the advantages of both strategies. The ability of high−order approximations is exploited to capture effectively the global response, while the superposed mesh provides an embedded approach suitable for detecting local effects.

Example 3: Snap−back of Cylindrical Panel
The static nonlinear response of a shell is discussed in this closing example. The test case is taken from [53] and considers an isotropic cylindrical panel loaded by a point force applied at its center and directed inward. The panel is illustrated in Figure 22, where the following parameters are considered: a = b = 50.8 mm, R = 254 mm, t = 0.635 mm. An isotropic material is considered with Young modulus E = 310.125 MPa and Poisson ratio ν = 0.3. The panel is under simply−supported boundary conditions along the straight edges, ξ 2 = 0, b, while free conditions are imposed along the curved ones, ξ 1 = 0, a.  Table 5. In all cases, the symmetry of the problem is exploited by modeling one quarter of the structure and imposing symmetry boundary conditions. A convergence study is performed by considering two levels of refinement for each model, as displayed in Figure 23. The h−refinement procedure is performed for a fixed polynomial order p = 1, with increased mesh resolution, as reported in Figure 23a. The p−refinement strategy is carried out by increasing the polynomial order, p = 1, 3, 6, on a mesh of 2 × 2 element, see Figure 23b. The ps−refinement is performed by increasing both the polynomial order p and the degree of s−refinement of the mesh s. More specifically, the p− and s−orders are increased simultaneously with steps of one (p = 1, 2, 3) and five (s = 0, 5, 10), respectively, see Figure 23c. For clarity, a summary of the FE refinement parameters is presented in Table 8 along with the number of degrees of freedom.  The results are presented in terms of load−deflection curves in Figure 24. An arc−length continuation procedure has been implemented for capturing the snap−back. As shown, all the solutions converge to the reference one provided sufficient refinement steps are performed. In particular, the results of Figure 24a demonstrate that 2 steps of refinement, leading to a total of 1327 dofs, are required for the h−model. On the contrary, one single step of refinement is needed for the p−model and the ps−L−model, as clear from the plots of Figure 24b,c. In these cases, the number of degrees of freedom is one order of magnitude lower for the p−model (343) and ps−L−model (252). These results demonstrate the advantage implied by the use of appropriate meshing and refinement strategies. Through a comparison of the p− and ps−refinement strategies with a classical h−refinement approach, the quality of the predictions is preserved, but the number of degrees of freedom can be as low as 1/5. This saving becomes even more important in the context of nonlinear analyses, where the overall process requires repeated matrix factorizations and linear systems to be solved.

Conclusions
This paper dealt with the application of the ps−version of the finite element method to the analysis of composite thin shells. The approach relies on the possibility of superposing multiple local meshes in arbitrary regions to refine the solution only where actually needed. No transition meshes are required nor special−purpose elements to connect areas with different grid densities.
The implementation proposed herein allows a wide class of structures to be analyzed, including shells with arbitrary curvature and innovative materials with variable in−plane elastic properties. Both geometrically linear and nonlinear features are introduced. Three different combinations of p− and s−refinements were considered, although other strategies could be easily developed and implemented.
The proposed advanced refinement technique is useful under several circumstances in the analysis of composite shells. The results illustrate this aspect for the case of local stress concentrations, where the possibility of superposing multiple layers of refinement can be exploited to achieve accurate stress predictions with reduced effort. For 2D laminated structures, advanced refinement can be of interest even in case of global responses, such as for the eigenanalysis of highly anisotropic plates. Indeed, owing to complex elastic couplings, the ability to capture local effects can have a drastic influence also on the predictions at a global level. For these problems, conventional methods based on high−order global approximations can be unsuitable if refined predictions are of concern. Among the different combinations of p− and s−refinements considered in the various test cases, no optimal choice can be established ex ante. The choice is, in general, problem−dependent.
The inherently hierarchical nature of the proposed tool suggests its use at different stages of the design process of composite shells. In the early design phases, simplified models with one or few macro elements with high−order interpolation can be the most suitable way for obtaining quick estimates. In more advanced phases, when increased detail is needed, the models can be easily improved by exploiting the local refinement capabilities to maximize the ratio between accuracy and number of degrees of freedom.