Composites Science and Technology 115 (2015) 43–51
Contents lists available at ScienceDirect
Composites Science and Technology journal homepage: www.elsevier.com/locate/compscitech
On the use of global–local kinematic coupling approaches for delamination growth simulation in stiffened composite panels R. Borrelli a, A. Riccio b, A. Sellitto b,⇑, F. Caputo b, T. Ludwig c a
Italian Aerospace Research Center (CIRA), Aerostructural Design and Analysis Lab., Via Maiorise, Capua, CE 81043, Italy Second University of Naples, Aerospace and Mechanics Engineering Department, Via Roma 29, Aversa, CE 81031, Italy c SMR SA, Blumenstrasse 14-16, CH-2502 Biel/Bienne, Switzerland b
a r t i c l e
i n f o
Article history: Received 2 February 2015 Received in revised form 8 April 2015 Accepted 11 April 2015 Available online 17 April 2015 Keywords: B. Delamination C. Finite element analysis (FEA) C. Multiscale modeling C. Stress transfer VCCT
a b s t r a c t This paper investigates kinematic coupling approaches for FE simulation of the mechanical behavior of complex composite structures. Two coupling methods belonging to the family of kinematic coupling approaches, namely the point-wise kinematic coupling and the weighted residual kinematic coupling are analyzed and assessed with respect to their application to non-linear structural FE analyses. The investigated kinematical coupling approaches have been implemented in the in B2000++ FEM code developed by SMRÒ and used to couple global–local domains characterized by different mesh refinements. A stiffened composite panel with an artificial bay delamination under compressive load has been adopted as numerical benchmark. Results in terms of stress distributions at domains interfaces, obtained with several FE models characterized by different coupling approaches and different mesh refinements, have been compared with the results obtained adopting a reference test-case model with mesh continuity and merged modes. Furthermore, the influence of the coupling approach on the delamination growth simulation has been assessed. Finally considerations about the computational effectiveness of the coupling approaches with respect to the reference test-case are introduced. Ó 2015 Elsevier Ltd. All rights reserved.
1. Introduction The use of composite materials is continuously increasing for aerospace applications due to their promising specific strength and stiffness. However, the complexity of composites damage phenomenology involving different failure modes, interacting each other, is seriously limiting the effectiveness of designing with these outstanding materials. One of the most common and critical failure modes is the delamination, which can have relevant effects on the reduction of the residual global strength of a composite laminate. Indeed, a delaminated composite panel under compression may experience delamination growth induced by local and global instabilities leading to the anticipated compressive failure. Finite element procedures for the simulation of the delamination growth induced by local and global instabilities are generally based on fracture mechanics concepts and involve the evaluation ⇑ Corresponding author. E-mail addresses:
[email protected] (R. Borrelli),
[email protected] (A. Riccio),
[email protected] (A. Sellitto),
[email protected] (F. Caputo),
[email protected] (T. Ludwig). http://dx.doi.org/10.1016/j.compscitech.2015.04.010 0266-3538/Ó 2015 Elsevier Ltd. All rights reserved.
of the Strain Energy Release Rates (ERR) at the delamination front [1]. Some numerical procedures use the total value of the ERR [2] without accounting for fracture modes separation while other techniques such as the Virtual Crack Closure Technique (VCCT) [3] allow to distinguish between the three basic fracture modes ERR contributions (GI, GII and GIII). The growth initiation is, generally, identified based on a growth criterion [4] (basically derived by fracture tests data and by the fracture toughness of the material GIC, GIIC and GIIIC). When the criterion is met along the delamination front, the delamination front is modified thus simulating the delamination opening. An example of linear delamination growth criterion presented in Eq. (1) [5]:
GI GII GIII þ þ ¼1 GIC GIIC GIIIC
ð1Þ
The ERR is a function of both the forces at the crack tip and the Crack Opening Displacements (COD). Hence, numerical tools for the prediction of delamination onset and propagation require a very fine three-dimensional FE model to accurately predict these quantities. On the other hand, a fine 3D discretization of complex structures is sometimes not feasible due to the increase of
44
R. Borrelli et al. / Composites Science and Technology 115 (2015) 43–51
calculation time. Therefore, with the aim of extending the field of applicability of such numerical tools to geometrically complex structures, global/local approaches are introduced [6–8]. In a global/local approach a very large number of elements is employed only in some critical areas of the investigated domain, also called local domains, where the solution is expected to show significant variations in small spatial and/or temporal intervals, while a coarser mesh is employed in the rest of the structure, also called global domain. The most common solution is to use transition elements to connect the refined local domain to the coarser global domain. An interesting example of transition element was developed by Davila and Johnson [9] to predict the compression strength of dropped-ply graphite–epoxy laminated plates. Each transition element was used to connect a shell element to a stack of brick elements in order to capture the fully three-dimensional response in the vicinity of the dropped ply. Numerical results were in good agreement with the experimental measurements. Cho and Kim [10] investigated the bifurcation buckling problem of a delaminated composite laminated by means of transition elements. In order to save computer resources without losing accuracy, layerwise elements were used in the local delaminated region. These elements were able to accurately describe the geometric deformations of the delaminated zone, being still computationally less expensive than three dimensional elements. In the rest of the laminate, where the deformation shapes are more regular if compared to the delaminated region, the first order shear elements were used. The transition element was then used to connect the global zone and the local delaminated one. The predicted buckling load was found in good agreement with the one obtained with a fully layerwise model. In [11] this transition element was demonstrated to be effective for the simulation of the post-buckling behavior of delaminated composites under compressive loads. However, the use of transition elements implies a relevant expertise in the formulation of the transition mesh because the presence of distorted elements can invalidate the solution in the transition regions. In order to effectively study local phenomena, more affordable tools can be used in FE codes to couple finite element models with different mesh densities and/or different element types. Such coupling methods allow carrying out very efficient and reliable global/local analyses on large structural components [12]. Alesi et al. [13] studied the behavior of skin/stringer interfaces in stiffened panels by a Multipoint Constraint (MPC) based method implemented in the ABAQUSÒ code [14]. Since the stress state at the skin/stringer interface is fully three-dimensional, brick solid elements were employed to model this critical local area; while a coarser two-dimensional mesh was used in the rest of the panel. The two different modeled areas were connected together via the MPC elements implemented in ABAQUSÒ. Krueger and O’Brien [15] and Krueger and Minguet [16] investigated the delamination propagation and the skin/stringer debonding respectively by shell-to-solid coupling elements implemented in ABAQUSÒ. Pietropaoli and Riccio [17] used the multipoint constraint approach implemented in ANSYSÒ code [18] in order to couple a coarser shell domain of a stiffened composite panel to a more detailed brick model representing a delaminated domain to test the effectiveness of a new tool for the inter-laminar and intra-laminar damage onset and evolution. The purpose of this paper is to investigate the behavior of a delaminated stiffened composite panel under compressive load and to test the effectiveness of two coupling methods implemented in B2000++ code developed by SMRÒ. The coupling methods belong to the kinematic coupling approaches family and are: The point-wise kinematic coupling.
The weighted residual kinematic coupling. In order to fully understand the peculiarity of these kinematic coupling approaches. Let’s consider two subdomains X1 and X2 connected with a common interface C, as shown in Fig. 1. The deformation field is indicated as x(X), where X are the coordinates of the reference configuration. The subdomains X1 and X2 are discretized such that the deformation fields, xe1 and xe2 , are C0-continuous across boundaries. The residual r(X) on X in C is equal to:
rðXÞ ¼ xe1 ðXÞ xe2 ðXÞ ¼ ue1 ðXÞ ue2 ðXÞ
ð2Þ
The conditions to be imposed on Eq. (2) result in a different coupling approach: In the point-wise coupling the residual is set to 0. In the weighted residual coupling the residual is minimized: r is integrated over C with a weight function and the minimum of the integral is obtained by setting the derivatives to 0 (zero). The couplings have been tested over a 3-dimensional finite element model representative of a stiffened composite panel with a defect located in the central bay and subjected to a compressive load; the couplings have been employed to connect the region of the defect (local domain), discretized with 3D elements, and the rest of the panel (global domain), discretized with shell elements. In Section 2 the test case is described while in Section 3 the two different kinematic coupling approaches are compared. 2. Test case and FE model description The numerical tests have been performed on the carbon-fiber reinforced panel, stiffened by three I-section stringers [19,20], schematically introduced in Fig. 2. An artificial circular delamination was located at the center of the bay and placed between the fifth and the sixth ply along the thickness of the skin. The stacking sequence of the skin was (+45°/45°/0°/90°)4s, while the stringers were made of 4 laminates stacked according to the sequence (4 5°/+45°/0°)2s. According to Fig. 2, the material 0° direction, which is also the loading direction, is parallel to the y axis. The mechanical properties of the adopted material system (T800/924) are summarized in Table 1 [21]. The experimental test has been performed under displacement control: the edge BC of the panel (see Fig. 2) has been clamped, while on the opposite edge (OA) a uniform displacement along the y-axis has been applied. The rotation on the edge OA has been constrained and the lateral edges (OC and AB) have been left ‘‘free’’. The area surrounding the delamination, which is of interest due to the high deformations and stress arising during the loading process, has been modeled by two layers of solid elements, while all the other regions have been modeled by shell elements. Global–local elements have been used to couple the two differently modeled regions. Contact elements have been also placed on the initial
Ω1
Γ
Ω2
Fig. 1. Subdomains X1 and X2 connected with the common interface C.
R. Borrelli et al. / Composites Science and Technology 115 (2015) 43–51
45
Fig. 2. Geometrical description of the numerical test case (units in mm).
Table 1 Material properties of the unidirectional composite lamina (T800/924). Property
Average (Cv)
Longitudinal Young’s modulus, E11 (GPa) Transverse Young’s modulus, E22 (GPa) Normal Young’s modulus, E33 (GPa) In-plane shear modulus, G12 (GPa) Out-of-plane shear modulus, G13 (GPa) Out-of-plane shear modulus, G23 (GPa) Poisson’s ratio, t12 Poisson’s ratio, t13 Poisson’s ratio, t23 Critical energy release rate mode I, GIC (J/m2) Critical energy release rate mode II, GIIC (J/m2) Ply thickness (mm)
155.21 (7%) 8.57 (3%) 8.57 (3%) 7.40 (5%) 7.40 (5%) 7.40 (5%) 0.36 0.36 0.36 280 575 0.125
delaminated area in order to avoid overlapping between the two delaminated sub-laminates, which could occur during the load application. A node to node penalty formulation has been adopted to reduce the computational effort. As results of a sensitivity analysis a penalty factor of 1 103 has been used, which has been found able to provide negligible penetration between contacting surfaces. Stringers have been modeled with shell elements connected to the panel via Multi-points Constraints (MPC).
Fail release fracture elements based on the Modified Virtual Crack Closure Technique (MVCCT [3]) and on the delamination growth criterion introduced in Eq. (1) have been implemented in the model, allowing the simulation of the delamination growth phenomenon. A shell/3D global local modeling approach has been used. Indeed a local area (80 mm 80 mm) including the delamination, which represents the area of interest where a detailed stress distribution evaluation is required, has been modeled with three-dimensional solid finite elements. For the rest of the panel a 2D-shell model has been employed. By means of such an approach, the accuracy of the three-dimensional model has been combined with the computational efficiency of a shell finite element model. In order to couple the shell domain to the solid one, shell-to-solid coupling elements based on kinematic constraints have been used. In particular, in the present paper kinematic coupling approaches have been investigated and compared: The point-wise kinematic coupling, implemented via Shell to Solid Coupling (SSC) elements which can be automatically created in the FEM research code B2000++ by using the ‘‘add_ss c_elements’’ directive.
Fig. 3. (a) Schematization of the solid finite element model in the delaminated area; and (b) Field Transfer boundary condition to the shell/solid interface.
46
R. Borrelli et al. / Composites Science and Technology 115 (2015) 43–51
Fig. 4. Schematization of investigated FE models.
Table 2 Summary of investigated FE models. FE model
Local mesh density
Kinematic coupling method
1 2 3 4 5 6
LMD1 LMD1 LMD1 LMD2 LMD2 LMD2
Point-wise Weighted residual n.a. – domains merged Point-wise Weighted residual n.a. – domains merged
The weighted residual kinematic coupling, implemented via the Common Refined Mesh (CRM) method and activated in B2000++ by the ‘‘add_field_transfer_coupling’’ directive. The solid (local) domain has been divided into three different zones. According to Fig. 3(a), where a quarter of the solid finite element model is considered, the three zones can be defined as: Zone I: is the initial delaminated area; contact elements are needed between the two sub-laminates in order to avoid overlapping. Zone II: is the zone where the delamination is allowed to grow. Fracture elements are introduced at the interface between the upper and the lower sub-laminates connecting corresponding pair of nodes. The connections between the nodes are released once the linear delamination growth criterion (Eq. (1)) is fulfilled. Zone III: is the zone where the solid mesh is connected to the shell mesh. The nodes at the interface between the two sub-laminates were ‘‘merged’’; hence the delamination has not been allowed to grow in this zone.
3. Comparative study: point-wise vs. weighted residual kinematic coupling approach In order to overcome some of the shortcomings (like stress jumps due to over-constraining), which involve the point-wise
Fig. 5. Comparison between numerical and experimental results: (a) LMD1 models; (b) LMD2 models; and (c) out-of-plane displacement computed at the center of delamination vs. applied strain.
R. Borrelli et al. / Composites Science and Technology 115 (2015) 43–51
47
Fig. 6. Delamination growth progression: (a) LMD1 models; (b) LMD2 models; and (c) delaminated area vs. applied strain.
coupling methods (such as B2000++ SSC elements), the Common Mesh Refinement (CRM) method has been developed and implemented in B2000++. This method belongs to the family of weighted residual kinematic coupling methods and it can be activated in B2000++ by using the ‘‘add_field_transfer_coupling’’ directive. The faces of solid elements and the edges of the shell elements are automatically coupled by the ‘‘field_transfer’’ boundary condition, which is automatically imposed by setting the ‘‘add_fie ld_transfer_coupling’’ directive. The surfaces where shell and solid elements are adjacent are automatically identified (see Fig. 3(b)). The point-wise coupling method and the weighted residual kinematic coupling method have been both applied to the investigated stiffened panel for comparison purposes. Basically, as
sketched in Fig. 4 and in Table 2, six finite element models have been introduced. Indeed, two levels of mesh density for the local domain have been considered. In particular, for the models identified as 1, 2 and 3 the lower mesh density (LMD1) has been used while for the models identified as 4, 5 and 6 the higher mesh density (LMD2) has been adopted. Moreover, three coupling approaches were investigated as described below: FE models 1 and 4: the local and the global domains have been coupled via SSC elements (point-wise kinematic coupling). FE models 2 and 5: the local and the global domains have been coupled via ‘‘Field Transfer’’ boundary condition (weighted residual kinematic coupling).
48
R. Borrelli et al. / Composites Science and Technology 115 (2015) 43–51
traditional approach for which no global/local elements are needed. It is worth to notice that in this traditional approach mesh refinements are needed also in the global domain in order to obtain matching meshes at the interface.
Fig. 7. Computational time comparison.
FE models 3 and 6: the local and global domains have been merged. For these last two models to allow the merge of the domains, the whole skin has been modeled with solid elements in order to have matching meshes at the interface between the global and local domain. This approach can be considered as the
In all the introduced FE models, fracture elements were used to simulate the delamination growth. The reaction vs. applied strain curves obtained with the LMD1 and LMD2 models are plotted in Fig. 5(a) and (b), respectively. In these figures, the curves obtained by using the point-wise kinematic coupling method, the weighted residual one and the traditional approach are compared each other and with the experimental results. No relevant difference can be appreciated between the different FE models and the global stiffness value of the panel correlates quite well the experimental one. In Fig. 5(c) the out-of-plane displacement of the lower sublaminate computed at the center of the delaminated area is plotted against the applied strain. The three numerical curves obtained with the LMD2 models are compared to experimental LVDT readings. Actually, the experimental test was performed in three steps in order to allow inspection of the defect at selected intervals. By analyzing the resulting three experimental curves, it seems that
Fig. 8. Out-of-plane displacement contour plot: (a) 2450 le; and (b) 3990 le.
R. Borrelli et al. / Composites Science and Technology 115 (2015) 43–51
the test was affected by LVDT’s calibration problems, hence, LVDT data have been introduced only for qualitative comparative purposes. The numerical curves obtained with the global/local models have been found very close to the one obtained with the traditional approach. The evolution of the artificial circular delamination for the LMD1 and LMD2 models is introduced in Fig. 6(a) and (b), respectively. Red fracture elements show the areas where the delamination growth criterion (Eq. (1)) is fulfilled and the separation between the initially bonded surfaces has been performed. No relevant difference has been appreciated, in terms
49
of delamination growth, for the investigated six different FE models. As a matter of facts, the weighted residual method is able to provide a more refined stress distribution at the global local interface with respect to the point-wise kinematic coupling approach. However, even if the stress distribution plays, in general, a critical role in damage evolution, inter-laminar progressive failure events, such as delamination growth events (which are the main focus of this paper) have been found almost insensitive to stress jumps at the global–local interface. This is the reason why, in terms of inter-laminar damage evolution no appreciable difference can be pointed out between the two analyzed approaches.
Fig. 9. (a) Stress (Syy) distribution progression for the LMD2 merged model; stress distribution at the interface: (b) 2450 le; and (c) 3990 le.
50
R. Borrelli et al. / Composites Science and Technology 115 (2015) 43–51
Indeed, influence of coupling approach on damage propagation is expected to be relevant for intra-laminar damages which have not been analyzed here. The delamination initiation phenomenon has been observed for an applied strain value of 3200 le. A relatively small angle (about 8°) has been found between the delamination propagation direction and the x direction, which means that the delamination propagates in a direction almost perpendicular to the load direction. For an applied strain of about 4266 le, the delamination reaches the boundary of zone II. Beyond this point the growth stops and the results of the analysis are no more reliable. This is a limit of these models without re-meshing, which are very effective from a computational point of view but do not allow the delamination growth region to change size during the analysis. The delamination evolution, in terms of shape, does not change significantly passing from LMD1 models to LMD2 models. However, the LMD2 models are able to provide a better resolution of the delamination shape as shown in Fig. 6(b). As already remarked, the kinematic coupling approach does not seem to affect significantly the results in terms of delamination growth. Moreover, the results obtained with global/local approaches are in very good agreement with the ones obtained with the ‘‘merged’’ model. In Fig. 6(c) the delaminated area is plotted against the applied strains for the following models: LMD1 model coupled via the ‘‘Field Transfer’’ boundary condition. LMD2 model coupled via the ‘‘Field Transfer’’ boundary condition. LMD2 model with the ‘‘merged domain’’. The delamination initiation load predicted by the three models is the same. Nevertheless, in the models with a more refined local mesh (LMD2), a slightly faster delamination grow is appreciated. The mesh dependency of the VCCT technique is a well-known trouble, which may be solved by tuning the mesh parameters with respect to an experimental DCB test, or applying mesh and time step independent approaches [22–26] which is not the scope of the investigated application. In Fig. 7 a comparison between the computational time needed for the investigated models is shown. The global/local approach allows to reduce the CPU time of about 50% (with respect to the corresponding merged model) without affecting the accuracy of obtained results. Models LMD2 coupled via SSC elements or via the Field Transfer boundary conditions seem to be a good compromise between reduced computational time and accurate results. In Fig. 8(a) and (b) the out-of-plane contour plots, predicted by the LMD2 merged model for the applied strains values of 2450 le and 3990 le, respectively, are reported. At the first strain value (Fig. 8(a)), the local buckling has occurred but the delamination defect has not yet started. At the second strain value (Fig. 8(b)), the delamination has instead significantly propagated. The progression of the Syy stress component along the section MN of Fig. 8(a) is plotted in Fig. 9(a). The stress distribution has been extracted at the mid-plane of the 3rd ply which is a 0° ply belonging to the upper sublaminate. It can be seen that, when the upper sublaminate starts to separate from the lower one, the load sustained by the buckled area starts to decrease. In this case, since the local domain was merged to the global domain, the continuity of the stress is guaranteed at the interface (x = 200 mm and x = 280 mm). In order to verify the performances of the two investigated coupling methods in terms of stress field continuity at the interfaces between the global and the local domain, the stress distributions computed by the model LMD2 coupled via SSC elements and the model LMD2 coupled via the Field Transfer boundary condition, have been compared to the one plotted in Fig. 9(a), which can be
Fig. 10. Summary of stress jumps at the interface between global and local domains.
assumed as the reference stress distribution, since it is not altered by the introduction of a global/local coupling approach. Details of the stress distribution comparison for the strain values of 2450 le and 3990 le are reported in Fig. 9(b) and (c), respectively. Stress jumps are obtained at the interface between the global and local domains (x = 200 mm) when using both the point-wise and the weighted residual kinematic coupling. Nevertheless, the stress jumps obtained with the weighted residual method are less relevant than the ones obtained with the point-wise method. Moreover, the extension of the area affected by ‘‘artificial’’ stress distribution and the deviation from the reference stress distribution are reduced too. Such smaller deviations introduces smaller errors when progressive failure analysis is activated. Stress jumps at the interface were computed also at a strain value of 1990 le (before local buckling) and for the Sxx and Szz stress component too. Such stress distributions are not reported for the sake of brevity but a summary of the stress jumps obtained with the point-wise method and with the weighted residual method is reported in Fig. 10. The weighted residual method allows reducing up to 87% the stress jumps at the interface with respect to the point-wise method. 4. Conclusions This paper shows that the results obtained with global/local approaches (point-wise or weighted residual kinematic coupling) correlate very well with those obtained with a traditional ‘‘domains merged’’ approach. Furthermore, the global/local models allows saving up to 50% of the computational time with respect to the traditional approach without losing in terms of results’ accuracy. It can be noticed that when using a kinematic coupling, the continuity of the displacement field is assured but stress jumps at the interfaces between the coupled domains are observed. The stress jumps obtained with the weighted residual method are less relevant (up to 87%) than the ones obtained with the point-wise kinematic coupling method. Moreover, the extension of the area affected by ‘‘artificial’’ stress distribution and the deviation from the reference stress distributions are reduced too. Such smaller deviations obtained with the weighted residual method introduces smaller errors when progressive failure analysis are performed. References [1] Hellan K. Introduction to fracture mechanics. Singapore: International student edition McGraw Hill; 1985. [2] Sun X, Tong L, Chen H. Progressive failure analysis of laminated plates with delamination. J Reinf Plast Compos 2001;20:1370–89.
R. Borrelli et al. / Composites Science and Technology 115 (2015) 43–51 [3] Krueger R. The virtual crack closure technique: history, approach and applications. Appl Mech Rev 2004;57(1–6):109–43. [4] Reeder JR. An evaluation of mixed-mode delamination failure criteria. NASATM-104210. [5] Reeder JR. 3D mixed-mode delamination fracture criteria – an experimentalist’s perspective. In: Proceedings of the 21st annual technical conference of the american society for composites, vol. 1, Lancaster, PA, 2006; p. 129–45. [6] Sellitto A, Borrelli R, Caputo F, Riccio A, Scaramuzzino F. Application to plate components of a kinematic global–local approach for non-matching finite element meshes. Int J Struct Integrity 2012;3(3):260–73. [7] Sellitto A, Borrelli R, Caputo F, Riccio A, Scaramuzzino F. Application of the mesh superposition technique to the study of delaminations in composites thin plates. Key Eng Mater 2012;525–526:533–6. [8] Sellitto A, Borrelli R, Caputo F, Ludwig T, Riccio A, Scaramuzzino F. 3D global– local analysis using mesh superposition method. Key Eng Mater 2013;577– 578:505–8. http://dx.doi.org/10.4028/www.scientific.net/KEM.577-578.505. [9] Davila CG, Johnson ER. Analysis of delamination initiation in postbuckled dropped-ply laminates. AIAA J 1993;31(4):721–7. [10] Cho M, Kim J. Bifurcation buckling analysis of delaminated composites using global–local approach. AIAA J 1997;35(10):1673–6. [11] Kim J, Cho M. Postbuckling of delaminated composites under compressive loads using global–local approach. AIAA J 1999;37(6):774–8. [12] Sellitto A, Borrelli R, Caputo F, Riccio A, Scaramuzzino F. Methodological approaches for kinematic coupling of non-matching finite element meshes. Procedia Eng 2011;10:421–6. [13] Alesi H, Nguyen VM, Mileshkin N. Global/local postbuckling failure analysis of composite stringer/skin panels. AIAA J 1998;36(9):1699–705. [14] ABAQUS Analysis User’s Manual – Part VI: Elements, p. 22.1.1–14. [15] Krueger R, O’Brien TK. A shell/3D modelling technique for the analysis of delaminated composite laminates. Composites: Part A 2001;32:25–44. [16] Krueger R, Minguet P. Skin-stiffener debond prediction based on computational fracture analysis. NASA/CR-2005-213915.
51
[17] Pietropaoli E, Riccio A. A global/local finite element approach for predicting interlaminar and intralaminar damage evolution in composite stiffened panels under compressive load. Appl Compos Mater 2010;18(2):113–25. [18] ANSYSÒ Mechanical, Release 14.0, Help System, Contact Technology Guide, ANSYS, Inc. [19] Greenhalgh E, Singh S, Nilsson KF. Mechanisms and modelling of delamination growth and failure of carbon-fibre reinforced skin-stringer panels. Composite structures: theory and practice. Philadelphia, PA: ASTM STP 1383; 2000. p. 49– 71. [20] Riccio A, Giordano M, Zarrelli M. A linear numerical approach to simulate the delamination growth initiation in stiffened composite panels. J Compos Mater 2010;44(15):1841–66. [21] Greenhalgh E. Characterization of mixed-mode delamination growth in carbon-fibre composites. PhD Thesis, Imperial College, London, 1998. [22] Pietropaoli E, Riccio A. On the robustness of finite element procedures based on virtual crack closure technique and fail release approach for delamination growth phenomena. Definition and assessment of a novel methodology. Compos Sci Technol 2010;70:1288–300. [23] Pietropaoli E, Riccio A. Formulation and assessment of an enhanced finite element procedure for the analysis of delamination growth phenomena in composite structures. Compos Sci Technol 2011;71:836–46. [24] Riccio A, Raimondo A, Scaramuzzino F. A study on skin delaminations growth in stiffened composite panels by a novel numerical approach. Appl Compos Mater 2013;20(4):465–88. [25] Riccio A, Raimondo A, Di Caprio F, Scaramuzzino F. Delaminations buckling and growth phenomena in stiffened composite panels under compression. Part II: a numerical study. J Compos Mater 2013;48(23):2857–70. http://dx.doi.org/ 10.1177/0021998313502742. [26] Riccio A, Raimondo A, Scaramuzzino F. A robust numerical approach for the simulation of skin–stringer debonding growth in stiffened composite panels under compression. Composites: Part B 2015;71:131–42.