Best Theory Diagrams for Laminated Composite Shells Based on Failure Indexes

This paper presents a novel approach to developing 2D structural theories for composite shells. The proposed approach uses the capabilities of the Carrera Unified Formulation (CUF) in conjunction with the Axiomatic/Asymptotic Method (AAM) to obtain the best theories for given structural layouts. Different structural cases are considered to analyze the influence of factors such as boundary conditions, lamination, and thickness on the accuracy of a model. The parameter chosen to evaluate a model’s performance is based on the Failure Indexes (FI) defined by the Hashin Failure criteria for unidirectional fiber composites. The outcome of this procedure is the Best Theory Diagram (BTD), containing the graphical representation of the highest accuracy as a function of the number of adopted unknowns. The results show the importance of higher-order terms to capture stress fields and the influence of thickness on the definition of the best theories.


Introduction
The development of advanced materials and their spreading adoption introduced new challenges in structural analysis concerning the accuracy and computational cost. The Finite Element Method (FEM) proved to be an extremely flexible and reliable tool for structural design. Different degrees of approximation can be chosen, opting for 3D elements or 2D/1D descriptions. In the last two cases, less accurate solutions must be expected, although at a significantly lower computational cost. Furthermore, reduced models can be enhanced by including a more complex description of the mechanical behaviors taking place over cross-section (1D) and the thickness (2D), albeit remaining efficiently solvable. This paper will focus on 2D formulations of composite shells.
Composite laminates are becoming increasingly predominant and used extensively in the aerospace and automotive areas. Starting from the Classical Lamination Theory (CLT), built from [1][2][3], and its later modification into the First-Order Shear Deformation Theory (FSDT) [4][5][6], many models were developed to obtain a more accurate description of the transversal behavior. This objective was achieved by increasing the number of generalized unknown variables, i.e., the nodal degree of freedoms (DOF), through higher-order polynomial thickness expansions [7][8][9][10] or the inclusion of non-polynomial terms [11][12][13][14]. The structural theories originated from this approach, often referred to as higher-order theories (HOT), approximate the laminate with an Equivalent Single Layer (ESL), meaning that the amount of adopted variables is independent of the number of layers. Furthermore, additional efforts were made to tackle numerical issues like the shear locking phenomenon [15][16][17] and overcome limitations relative to the interlaminar continuity of displacement and transverse stresses (C 0 z -requirements). Zig-zag [18] and Layer-Wise (LW) [19][20][21][22] models were thus developed, together with the new mixed finite element formulations [23,24], like those based on the Reissner Mixed Variational Theorem (RMVT) [25]. Another class of models stemmed from the Asymptotic approach [26] with the advantage of the a-priori determination of the accuracy with respect to 3D solutions.
This paper obtains higher-order theories from the Axiomatic/Asymptotic Method (AAM) [27][28][29]. This process & M. Petrolo marco.petrolo@polito.it starts by axiomatically choosing a maximum order of the polynomial expansions and assuming the full model as the reference. Terms are then suppressed, and for each expansion, a level of accuracy can be evaluated with respect to the reference solution. A fundamental aspect of this approach is the definition of the accuracy parameter. As shown in previous works [30][31][32], depending on the choice of the output, different theories can emerge as the optimal ones for a given structural problem. This paper considers a new parameter for selecting the best theories for 2D models of composite shells. It is based on the failure indexes evaluated according to the Hashin Criteria [33]. Different structural problems are considered to highlight the effect of features such as thickness, boundary conditions, curvature, and stacking sequence. Implementing the AAM is done in the framework of the Carrera Unified Formulation [18], an efficient and generalized approach to obtain the governing equations for virtually any finite element model. This paper is organized as follows: Sect. 2 presents CUF and FEM formulations; Sect. 3 introduces the Hashin Failure Indexes; Sect. 4 describes the AAM approach; numerical results are provided and discussed in Sect. 5; conclusions are drawn in Sect. 6.

Carrera Unified Formulation and Finite Elements
The Carrera Unified Formulation (CUF) [34] is used to efficiently derive all the theories required for the analyses presented. Using the reference system employed in Fig. 1, the displacement field can be expressed as The Einstein notation acts on s, while F s are the thickness expansion functions. u s is the vector of the generalized unknown displacements and M is the number of expansion terms. The following displacement equations describe a complete fourth-order model (E4) using Taylor-like polynomial expansion: The order and type of expansion are free parameters. This paper uses Equivalent Single Layer (ESL) formulations, meaning that a single set of variables is assumed across the whole thickness of the structure, independently from the number of actual layers. Throughout the paper, some comparisons will be performed with models adopting a Layer-Wise (LW) formulation. Concerning the geometrical relations, metric coefficients H k a , H k b and H k z of the k th layer are defined in the case of multilayered shells, As shown in Fig. 1, R k a and R k b are the principal radii of the middle surface of the k th layer, A k and B k the coefficients of  the first fundamental form of X k . In this paper, only constant radii of curvature were considered, A k = B k = 1. Strains can be written as where, For the stress-strain relations, it follows that: where: The FE formulation adopted in this paper uses a nine-node shell element (Q9) based on the Mixed Interpolation of Tensorial Component (MITC) method [17]. The displacement vector becomes: with u si and du sj being the nodal displacement vector and its virtual variation, respectively. By substituting them in the geometrical strain expressions described by Eq. 4, one obtains: MITC contrasts the membrane and shear locking via a specific interpolation strategy for the strain components on the nine-node shell element.
Strains aa m1 , bb m2 , c ab m3 , c az m1 , and c bz m2 stem from 10. Subscripts m1, m2, and m3 indicate quantities evaluated at the point groups (A1, B1, C1, D1, E1, F1), (A2, B2, C2, D2, E2, F2), and (P, Q, R, S), respectively, see Fig. 2. The interpolating functions are Lagrangian functions grouped in the following arrays: Table 2 Example of best theory representation DOF u a1 u b1 u z1 u a2 u b2 u z2 u a3 u b3 u z3 u a4 u b4 u z4 u a5 u b5 u z5 %E AVG Via the Principle of Virtual Displacements (PVD) for the static analysis, the equilibrium equation reads The 3Â3 matrix k k ssij is the fundamental mechanical nucleus whose expression is independent of the expansion order. p k sj is the load vector. For a more detailed description of the FE derivation process, the reader can refer to [35].

Hashin Failure Criteria
For the present study, the Hashin failure criteria [33] is adopted. This method, developed for unidirectional fiber composites, is based on a quadratic approximation of the failure-free envelope in the stress space. Its definition relies on using the transversely isotropic invariants of the stress tensor written in the material reference system. Hashin's approach identifies two main failure mechanisms, one tied to fibers and one to the matrix. For each, a tensile and a compressive mode can be distinguished. The mathematical formulation thus relies on the following four equations: Tensile Fiber Mode (r 11 [ 0): Compressive Fiber Mode (r 11 \0): Tensile Matrix Mode (r 22 ?r 33 [ 0): Table 5 Best theories, SS, a/ h=100 Compressive Matrix Mode (r 22 þ r 33 \0): X T;C and Y T ;C are the tensile and compressive strength for fibers and matrix, respectively, while S 12 , and S 23 represent the material shear strengths. Concerning the evaluation of the failure indexes, all the results presented here refer to the strength properties of the IM7/977-3 laminate, gathered from [36] and summarized in Table 1.

Axiomatic/Asymptotic Method and Best Theory Diagrams
This section presents a method to select higher-order generalized variables and build structural theories called the Axiomatic/Asymptotic Method (AAM). For the sake of brevity, only polynomial expansion terms are considered primary variables. The preliminary, axiomatic step of the AAM is the choice of the maximum order of the expansion to be reached. This work considered fourth-order models as they usually provide high accuracy. This choice results in a total of 2 15 possible combinations of active terms due to the fifteen available generalized unknowns that can be used. Of those fifteen terms, the three zeroth-order ones u a1 , u b1 , and u z1 were always kept, reducing the number of possible Table 6 Best theories, SS, a/ h=10 combinations from 2 15 to 2 12 . This choice was motivated by the considerable impact of these unknowns on the model's accuracy, making them essential for deriving meaningful results. Moreover, the full fourth-order expansion model (E4) is the reference to compute errors. The asymptotic aspect of the procedure requires the evaluation of the performance of each model, and this can be done in different ways, depending on the outputs of interest. As described in previous works [30][31][32]37], the performance parameter can be defined as a percentage error over a single displacement at a specific point, or as an average of the errors relative to a certain number of modal frequencies obtained through free vibration analyses. In this paper, the quality of a single theory is defined by the average error over the two failure indexes, namely, the Tensile Matrix Mode (f MT ) and the Compressive Matrix Mode (f MC ) indexes, evaluated at the top and bottom of the central point of the shell. The per-centage errors relative to the E4 model were evaluated as follows: The parameter chosen to assess the quality of a model is then evaluated as: The outcome of this selection procedure is provided using the Best Theory Diagram, as in Fig. 3.
The BTD provides the errors for different numbers of active variables obtained via models using specific sequences of expansion terms. The resulting Best Theories can be conveniently represented in tabular form for every set of structural parameters. The adopted graphical Table 7 Best theories, SS, a/ h=5 convention uses black and white triangles to indicate active and inactive expansion terms, respectively, as in Table 2.
For comparisons purposes, a second performance parameter was also considered, the percentage error over the maximum transverse displacement [38],

Numerical Results
Several structural cases of square doubly-curved shells were considered, all gathered from [38]. A total of seven different combinations were selected, varying in thickness, stacking sequence, and boundary conditions. Curvature radii were always kept equal to R/a=5. The material properties common to all the cases analyzed in this paper are reported in Table 3.
The load is a bi-sinusoidal pressure distribution acting on the upper surface, defined as p z ¼p z sinðpa=aÞsinðpb=aÞ, with unit amplitudep z =1. Unless otherwise specified, a quarter of the shell was used with a 4Â4 mesh of Q9 elements and symmetry boundary conditions. Similarly to [38], the results were obtained using an ESL formulation, and the reference is a full fourth-order E4 Taylor expansion.

Preliminary Assessment
The verification considered the axial stress and compared it to [39]. The non-dimensionalized stress r aa was used, The entire geometry discretized with a 9Â9 mesh was adopted for the stress evaluation, following the one used in [39], together with an ESL formulation and full fourth-order Taylor expansion. The reference paper contains both the E4 solution and the one obtained with a fourth-order Layer-Wise description, referred to as L4. Table 4 has the evaluated stresses, and the two reference solutions brought using the E4 and L4 formulations. A perfect match was found with the results provided in [38] for all considered configurations. For thinner shells, the E4 and L4 models provide identical results. However, the gap between the two solutions widens with an increase in the a/h parameter. This phenomenon is related to the necessity of a much more detailed description of the transverse behavior typically required for the analysis of thicker structures, here provided by the L4 formulation. However, the difference between L4 the E4 remains small, and E4 was kept as a reference for all the following cases.

Best Theory Diagrams, [0 /90 /0 ]
The first presented case group consists of the simply-supported (SS) spherical shells with symmetric lamination [0 / 90 /0 ]. Figures 4, 5, and 6 are the BTD for shells with thicknesses a/h=100, 10, and 5, respectively. Tables 5, 6, and 7 show their corresponding Best Theories, paired with the achieved average error. From them, different displacement field descriptions can be retrieved; for instance, the best model with ten DOF in the case of an SS shell with a/ h=10 is Furthermore, the distributions of the r aa and r az stresses along the thickness were analyzed, comparing the best models with 5, 10, and 15 active DOF. The aim is to assess the best models obtained via FI computed in a limited number of points in obtaining the through-the-thickness distribution of stresses. The non-dimensionalized form for r az is r az ¼ r az ð0; a=2Þ ĥ p z a ð24Þ The results from a third-order Layer-Wise model (L3) were also added. Figure 7 shows the r aa and r az for the a/h=100 case, Fig. 8 for a/h=10, and Fig. 9 for a/h=5. Table 8 compares the active terms for some of the Best Theories obtained for the three symmetric lamination cases, specifically for seven, nine, and eleven total degrees of freedom. The expansion terms are presented from lowest to highest order going from left to right in each sub-table row. The same graphical convention previously introduced to indicate an active/inactive term is used. Table 9 compares instead the Best Theories for the same number of active DOF obtained using the error over the central transverse displacement. These results were gathered from [38].
The following considerations stem from the results above, • Generally, the BTD obtained using the average error over the FI have relatively higher errors than those stemming from the error over the vertical displacement. The gap tends to be lower for thinner structures. Best models with six or fewer DOF have significant errors in the evaluation of FI due to the influence of higher-order generalized displacements over the stress distributions.
• The Best Theories found for the thinner configurations have zeroth-and first-order terms, with in-plane components -u a2 and u b2 -being predominant. These are followed by the second-order u z3 and the first-order u z2 . u b terms above the second-order appear to be the least influential for this case.
• For the moderately thicker configuration with a/h=10, first-order u a2 and u b2 remain the most relevant terms. u b terms of second and third-order are more frequent than the a/h=100 case, together with the fourth-order u z5 .
• Third-order terms are the most frequent for the thicker shell with a/h=5. More in detail, the linear term u a2 is almost wholly neglected, being replaced by the thirdorder u a4 . The effect of the inclusion of a cubic term was particularly beneficial for the description of the stress distribution along the thickness of the structure. As shown in Fig. 9a, this allowed overcoming the limitation of a constant stress distribution imposed using an FSDT model.
• On the effect of using a particular selection parameter, Tables 8 and 9 show that for small numbers of active DOF, both methodologies produce very similar, if not identical, best theories. The increase in active terms and thickness can make some noticeable changes between the results stemming from the two criteria. For thinner     configurations, the best models mainly differ from second or higher-order terms, while, in the thicker cases, differences can already be found in active first-order terms. Some of the most significant variations are the loss of the second, third, and fourth-order u b terms for the a/h=100 case compared to the best theories obtained with the error over the central transverse displacement Table 12 Best theories, CF, a/ h=100, [0 /90 /0 /90 ] DOF u a1 u b1 u z1 u a2 u b2 u z2 u a3 u b3 u z3 u a4 u b4 u z4 u a5 u b5 u z5 %E AVG Table 13 Best theories, CF, a/ h=5, [0 /90 /0 /90 ] DOF u a1 u b1 u z1 u a2 u b2 u z2 u a3 u b3 u z3 u a4 u b4 u z4 u a5 u b5 u z5 %E AVG  and the greater relevance of fourth-order terms for the a/ h=5 case in theories evaluated through FI.
• As shown in previous works [30], there may be oscillations in the error/DOF plot, i.e., an increasing number of generalized unknown variables does not necessarily lead to higher accuracy. Adding single terms may be detrimental if all those with similar influence are not retained.
• Overall, using ten DOF leads to through-the-thickness distributions similar to the full E4 case. However, for transverse shear stress, as known, LW models provide more accurate results.

Best Theory Diagrams, [0 /90 /0 /90 ]
Four shell configurations with asymmetric lamination [0 / 90 /0 /90 ] are presented in this section. These consist of two SS cases with a/h=100 and 10, followed by two clamped-free (CF) ones with a/h=100 and 5; the clamped edges are those parallel to the b-direction. Tables 10 and 11 show the Best Theories obtained for the SS shells. Tables 12  and 13 are those for the CF cases. For instance, the best model with ten DOF for a CF shell with a/h=5 is   For these two latter cases, the plots of r aa and r az along the thickness are presented in Figs. 16 and 17. The comparison between several best theories with the same number of active DOF across the four considered configurations is shown in Table 14, while  Table 15 presents the best models obtained in [38] through the transverse central displacement.
From the second set of results, it can be observed that • Starting from the SS cases, linear terms u a2 and u b2 are still the most relevant. For the thinner configuration with a/ h=100, second-order terms are more frequent than thirdorder ones, while the opposite is true for the thicker case with a/h=10. By comparing the results for the asymmetric lamination with their symmetric counterparts, it can be noted that there is a slight reduction in the relevance of thirdorder terms for the thinner shell in favor of some fourth-order ones. Fourth-order in-plane terms gain more weight for the thicker case, together with the third-order u b4 one.
• For the CF configurations, third-and fourth-order terms become more dominant, especially for the in-plane terms u a4 and u a5 . For the a/h=5 case, linear terms barely present up to nine active DOF, replaced by third-order in-plane ones. These expansion terms allow even simpler models to provide solutions much closer to the reference one, as evident from the transverse stress distribution presented in Fig. 17a.

•
The comparison provided in Table 14 summarizes the highlighted behaviors: high-order expansion terms are mandatory in the case of thicker structures.

•
The best theories derived using FI present some differences from those obtained through the transverse displacement. Starting from a small number of active DOF, those relative to SS cases are identical. At the same time, more significant variations can be found for the clamped-free cases. By increasing the number of expansion terms, the deviations are few and mainly pertain to third-or fourth-order terms.

Conclusions
This paper presented a novel approach to building shell models for the static analysis of composite structures. Failure indexes were used to evaluate the accuracy of reduced models and investigate the influence of higher-order generalized displacement variables via the Axiomatic/Asymptotic approach. Equivalent Single Layer models were considered with up to fourth-order terms. The Carrera Unified Formulation was used to obtain the governing equations and related finite element arrays. Different combinations of boundary conditions, stacking sequences, and thickness values were considered. Best Theory Diagrams were obtained in which, for a given number of nodal DOF, the set of generalized displacement variables to use to maximize the accuracy can be read. The following conclusions can be drawn: • Higher-order terms are essential to determine Failure Indexes and mandatory for thicker shells, where the role of third-order terms is decisive.
• The proper choice of higher-order terms can lead to models requiring very few DOF, e.g., for thicker shells the distributions of shear along the thickness do not present significant variations as soon as third-order terms are considered. However, the sole use of Equivalent Single Layer higher-order models is insufficient for high-fidelity prediction of transverse stresses.

•
The effect of boundary conditions and stacking sequence was negligible compared to the thickness, mainly affecting the activation of different expansion terms of the same order. For shells with thickness a/h=100 and 10, the first-order terms u a 2 and u b 2 are always present, with third and fourth-order ones being slightly less beneficial for the thinner cases.
The development of structural models is problem-dependent, and BTD may provide guidelines on those variables with more significant influence or indicate the accuracy of a given structural theory. However, the construction of BTD may be computationally expensive due to the strong dependence on input parameters. The computational burden may be avoided, and future works will explore the simultaneous use of multiple performance parameters to improve the robustness of the reduced models and use machine learning techniques to substitute FEM analyses and obtain the best models by trained networks.
Funding Open access funding provided by Politecnico di Torino within the CRUI-CARE Agreement.

Declarations
Conflict of interest No conflict of interests/competing interests to declare.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not Best Theory Diagrams for Laminated Composite Shells based on... 217 included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons. org/licenses/by/4.0/.