# Effect of delamination defects on buckling and growth characteristics in composite laminates

S.S. Venkat, S. Scheffler, P.M. Anilkumar, E. Baranger, R. Rolfes

► To cite this version:

S.S. Venkat, S. Scheffler, P.M. Anilkumar, E. Baranger, R. Rolfes. Effect of delamination defects on buckling and growth characteristics in composite laminates. Composite Structures, 2025, 369, pp.119245. 10.1016/j.compstruct.2025.119245. hal-05360315

HAL Id: hal-05360315

https://hal.science/hal-05360315v1

Submitted on 12 Nov 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Effect of delamination defects on buckling and growth characteristics in composite laminates

S.S. Venkatª,*, S. Schefflerª, P.M. Anilkumarª,*, E. Barangerᵇ, R. Rolfesª

<figure><img src="image_2.png" /></figure>

## ARTICLE INFO

### Keywords:

- Delamination defect

- Buckling

- Imperfection analysis

- Finite elements

## ABSTRACT

Delaminations pose a critical threat to the safety and reliability of composite structures due to increased risk of buckling and further growth under compressive loads. Understanding the interplay between how even seemingly minor interface defects can amplify the risk of buckling and how buckling in turn can affect the growth characteristics of existing delaminations is therefore crucial during the design phase of the structure. In this study, we focus on the column buckling case to get an understanding of growth characteristics of pre-existing delamination under compression. A quasi-static "Effect of Defect" study has been conducted using geometrically nonlinear finite element (FE) simulations. The study has been performed by maintaining a constant applied load while incrementally increasing the delamination length in the models. Different through thickness delamination configurations have been investigated in a generic quasi-isotropic laminate and the configuration with delamination near the laminate surface has been found to be most critical. It is observed that as the delamination length is increased, the laminate initially experiences local buckling. This is followed by unstable growth and a mixed local/global buckling failure when the delamination length reaches a critical size. For this near-surface configuration, the influence of laminate dimensions has been investigated and a comparison has been made between the original quasi-isotropic laminate and a stiffer, zero-dominant laminate. It has been demonstrated that in the stiffer laminate, under the same applied load, local buckling occurs without leading to unstable growth, even for very large delamination sizes. Therefore, in this case, if growth occurs under a service load equalling the applied load, it is expected to remain stable without transitioning to unstable growth. By understanding these distinct growth characteristics within the context of delamination buckling, designers can more effectively develop inspection strategies that improve the monitoring of damage growth.

## 1. Introduction

Fibre reinforced composites, owing to their high specific stiffness, strength and corrosion resistance, have found diverse applications, especially within the aerospace and wind energy sectors. These composites are usually produced by stacking plies of continuous fibre-reinforced plastic matrix in various orientations. The excellent in-plane properties imparted by the reinforcements are however not often matched by the through-thickness properties of laminates, making them prone to failure at ply interfaces known as delaminations. The troubling nature of delaminations stems from the difficulties in detecting them and the severe detriment they can cause to the mechanical performance of the structure [1]. As a result, they remain a major concern for designers. Delaminations can result from impact loads or be induced by both static or fatigue loads, due to factors such as matrix cracks, stress concentration near material or geometric discontinuities, or the

growth of interfacial defects introduced during manufacturing [2–4]. Under compression loads they can significantly decrease the stiffness and load-carrying capacity of composite structures, by potentially triggering premature buckling and subsequent delamination growth [2,5]. Therefore, it is crucial for the designer to address this issue during the design phase to meet the certification requirements. In this regard, Pascoe [6] highlighted two damage tolerance concepts for designing the composite structure. The currently more prevalent ‘no growth’ approach involves keeping in service loads below fatigue thresholds to prevent the growth of existing damage, such as delaminations. However, this method is conservative and therefore imposes a weight penalty. Contrary to the ‘no growth’ damage tolerance concept, the ‘slow growth’ concept allows existing damage originating from impacts or manufacturing imperfections in composite structures to grow under service loads such as fatigue. This is permissible as long as the growth

<footnote>ªInstitute of Structural Analysis, Leibniz Universität Hannover, Appelstraße 9A, 30167 Hannover, Germany</footnote>

<footnote>ᵇUniversité Paris-Saclay, CentraleSupélec, ENS Paris-Saclay, CNRS, LMPS - Laboratoire de Mécanique Paris-Saclay, 91190, Gif-sur-Yvette, France</footnote>

<footnote>* Corresponding author.</footnote>

<footnote>E-mail address: s.venkat@isd.uni-hannover.de (S.S. Venkat).</footnote>
remains slow, stable, and predictable, and the residual strength stays above the design limit load (DLL). The need for damage to be slow and stable is to define meaningful inspection window and it should predictable to enable designers to establish effective inspection schedule with enough checks within that window. A key aspect of applying the slow growth concept to delamination buckling is assessing whether the growth of an existing delamination could be unstable under a given load and laminate configuration, and if so, identifying the critical delamination size associated with it. To conduct such an extensive parametric study on growth characteristics, analysis methods are a requirement. The following paragraphs outline the methods used in the literature to study the effects of delamination on buckling and post-buckling behaviour, highlighting the advantages and limitations of each method. This comprehensive review is essential to identify the most suitable method for this study.

Although numerous studies have explored this topic, with early works focusing on the influence of delamination on the critical buckling load and mode shape, they did not account for the subsequent growth of delamination [7–16]. However, few of these studies address whether the growth is expected to occur, and if so, what its expected shape might be [7,11]. For instance, an early study by Whitcomb and Shivakumar [7] analysed the growth characteristics of a near-surface rectangular embedded delamination in a composite plate, caused by thin sublamine buckling. This was done by calculating the energy release rate (G) at the delamination front using the Virtual Crack Closure Technique (VCCT). It was observed that for a square and rectangular embedded delamination, there is a large variation of G along the delamination front due to which a self-similar growth is not expected. Additionally, the study demonstrated that possibility of delamination growth in the loading or transverse direction is influenced by factors such as the aspect ratio and absolute size of delamination. Similarly, Hwang and Huang [11] conducted a nonlinear finite element (FE) analysis to investigate buckling and assess the likelihood of postbuckling delamination growth in composite laminates with two through-width delaminations of varying lengths at different through-thickness locations. Multiple analyses were conducted by varying the length ratios of the two delaminations. By calculating G at the delamination front using VCCT, they were able to identify which of the two delaminations is likely to grow and also determine the associated fracture mode (opening, shearing or mixed mode).

Conversely, other studies have focused solely on buckling and post-buckling deformation behaviour using computational models, without discussing about the growth. For example an analytical study by Gu and Chattopadhyay [9] investigated the influence of single through-width delamination. From their experimental findings, it was reported that the through-thickness location of the delamination influences the mode of buckling. Their analytical model which comprised of higher order laminate theory provided a good agreement of buckling stress and postbuckling out-of-plane displacement path with experiments. In numerical and experimental works by Hwang and Liu [10] and Aslan and Şahin [13], influence of through-width multiple delaminations was studied in detail. The through-thickness distribution featured the longest delamination near the surface, tapering to the shortest delamination at the midplane, to replicate the effect of impact damage. Hwang and Liu [10] utilised a nonlinear FE based methodology to predict the postbuckling out-of-plane displacement growth, while Aslan and Şahin [13] employed linear bifurcation analysis (LBA) to predict the initial buckling event and its corresponding load. Both the studies demonstrated that the initial buckling event was solely controlled by the longest near-surface delamination, while the shorter underlying delaminations had no influence on it. Hwang and Liu [10] further illustrated that as the applied load was increased, multiple local buckling events corresponding to shorter delaminations beneath the longest near-surface delaminations occurred. Eventually, this transitioned into a mixed mode when the thickest

lowest sublamine buckled in opposite direction. Furthermore, in several studies, researchers [12,14,15] have developed analytical methods to predict the buckling and postbuckling deformation behaviour of delaminated laminates and compared these methods to nonlinear FE solutions. In their initial study [12], they employed classical laminated plate theory (CLPT) and developed a formulation based on the Rayleigh–Ritz technique to investigate single near-surface delamination. In their subsequent study [14], they compared CLPT with first-order and higher-order shear deformation theories (FSDT and HSDT) for the same problem. Their final study [15] focused on examining the effects of two through-width delaminations and a single circular embedded delamination using FSDT. Although the methods presented in the previously discussed literature offer initial insights into the problem of delamination buckling, they are insufficient for understanding the growth characteristics.

In general, there are wide range of numerical studies considering the delamination growth and they are well documented in literature [2,5,17–25]. For the analysis of postbuckling deformation, alongside nonlinear FE analysis [2,5,17,19–25], other approach like layerwise B-spline finite strip method (FSM) [18] has also gained attention over the years. In these methods the delamination growth has been incorporated using techniques such as VCCT with release of constraint between the tied surfaces [17–19] and cohesive zone model (CZM) [2,5,20–22,24].

Gaudenzi et al. [17] developed a modification to the incremental continuation method (by Riks technique) to improve the numerical FE analysis of unstable delamination growth in delaminated laminates. They analysed laminates with single through-width and circular embedded delamination. For the through-width delamination, a 2D plane strain model (with geometry and loading invariant along the width) was utilised, while a 3D model was employed for the embedded delamination. The delamination growth which was simulated using VCCT, was predicted to be significant for through-width configuration but only marginal for embedded configuration. Comparison with experiments in literature yielded a good agreement. In another study Wang and Zhang [18] employed B-spline FSM for deformation analysis and VCCT for delamination growth analysis of laminates with through-width delamination. By comparing local laminate surface strains with experiments conducted by Kutlu and Chang [26], they demonstrated that their method yields a good prediction of buckling loads corresponding to thin and thick sublamine. Similarly, Liu et al. [19] used VCCT approach to compute energy release rate for delamination growth in a geometric nonlinear FE analysis and made a comparison of their prediction with the results by Wang and Zhang [18]. Although VCCT has shown good agreements for these simple 3D cases and 2D plane strain problems, they still have difficulty propagating cracks when the crack front is not aligned with the mesh in 3D, and require additional techniques such as automated and robust remeshing procedures [27] or gradual release of forces [28]. Moreover, such linear elastic fracture mechanics (LEFM) techniques cannot predict crack nucleation and needs an initial crack in the model.

Alternatively, researchers have utilised CZM within a geometrically nonlinear finite element framework to simulate postbuckling delamination growth in delaminated laminates [2,5,20–22,24], as well as to analyse stiffened panels under compression [29]. In these studies fracture mode-dependent traction-separation law was employed for CZM. In one such studies, Mekonnen et al. [2] demonstrated that altering the through-thickness location of delamination from near-surface to midplane influences the mode of buckling, and consequently, affects the delamination growth load as well. It was also reported that after thin sublamine buckling, the sequence of events involving thick sublamine buckling and delamination growth can vary based on the delamination length and its through-thickness location (up to intermediate locations). In another work by Nikrad et al. [20], an analytical model based on layerwise approach within the framework of FSDT was developed and combined with VCCT for growth capability and the results were compared with the FE predictions using CZM.
The comparison of postbuckling out-of-plane displacements between the analytical and numerical models showed good agreement for edge and through-width delamination configurations. Additionally, it was demonstrated that off-centre through-width delamination can produce asymmetric instability configurations and may also trigger asymmetric delamination growth. Furthermore, CZM has also been employed to study the growth of embedded delamination under compression [21] and due to snap-buckling of sublaminate under flexural loading [21]. Although accurate energy dissipation in CZM requires multiple elements in the fracture process zone [30], CZM does not require an adaptive remeshing technique around the growth front and can predict both crack initiation and propagation.

The literature reviewed above provides an overview of methods used to model post-buckled delamination growth behaviour and also offers insights from parametric studies on the effects of delamination shape, size, and location on the force-displacement response. However, the focus of these studies has not been on analysing whether the delamination growth can be unstable for a given laminate configuration under an applied load. The growth is said to be unstable if there is significant extension of a pre-existing delamination under a given constant load. In order to study growth characteristics, a quasi-static ‘Effect of Defect’ study has been performed in this work on a laminate with through-width delamination. This study involved a series of analyses in which a constant load has been applied while the delamination length has been progressively increased in the models. The aim is to determine if the growth of an existing delamination can become unstable and if so, at what critical length. The research examines various through-thickness locations of delamination to identify the most critical configurations, and also explores whether unstable growth can occur when factors such as load level, laminate dimensions or laminate layup are altered. An essential part of implementing the ‘slow growth’ damage tolerance concept to delamination buckling is evaluating whether delamination growth can become unstable during the operational lifespan under specific service loads. By exploring the factors that trigger unstable delamination growth in laminates, this study therefore paves the way towards it. A nonlinear FE framework, coupled with CZM has been employed to analyse the evolution of out-of-plane displacement due to buckling and delamination growth.

The paper is organised in four sections. In Section 2 the methodology has been discussed which comprises the general aspects of the study, description of CZM, modelling framework and a validation example. The subsequent Section 3 presents and discusses the results of the parametric study. The paper ends with a section on conclusions along with the outlook in Section 4.

## 2. Methodology

### 2.1. General aspects of the proposed study

The FE modelling has been carried out using a commercially available software package, ABAQUS 2023. For an applied load of constant magnitude, the evolution of out-of-plane displacement with increasing delamination length and the critical delamination length corresponding to unstable growth have been examined. Various through-thickness locations of delamination have been studied in a generic example of a quasi-isotropic (QI) laminate with a layup of [0$_{3}$/ - 45$_{3}$/45$_{3}$/90$_{3}$]$_{S}$. This QI laminate has been selected for its simplicity, featuring just 8 sublaminates (eg: 0$_{3}$). Each sublaminate can be modelled with a single element through the thickness, helping to minimise calculation time.

Fig. 1 illustrates further details of the laminate and the investigated delaminated configurations. The configurations are chosen such that they are representative cases for composite delamination defects in literature. In a study by Donadon and Almeida [31], it is mentioned that delaminations typically occur between two layers of dissimilar orientations within a laminate. Another study on barely visible impact damage by Baluch et al. [32] indicates that largest delaminations are often

observed at interfaces farthest from the impact point, near the laminate surface on the opposite side, due to a phenomenon known as back-face splitting. These findings suggest that near surface and intermediate locations in this study are plausible locations for delaminations. Additionally, a mid-plane delamination at the interface between layers of the same orientation has been chosen, as it creates two sublaminates with identical stiffness, which may influence the buckling mode and, in turn, affect the critical delamination length associated with unstable growth. The laminate dimensions and other details in Fig. 1 have been selected based on findings from the literature [2,5] and preliminary investigations. The study also aims to provide a comparative analysis with laminates of different dimensions to assess the influence on growth characteristics. Finally the influence of laminate stiffness has also been analysed. To achieve this, the layup of the aforementioned QI laminate has been modified to increase stiffness while maintaining the same thickness. Additionally a brief discussion has been presented on the influence of the magnitude of the applied load in the stiffer laminate. An important step in this study is to simulate the delamination growth and it has been discussed further in the next subsection.

### 2.2. Static cohesive zone model

The buckling mode of delaminated laminates under compression can vary based on the through-thickness location and length of the delamination. The buckling mode in turn directly influences the fracture mode associated with the subsequent growth of delamination. Considering the same, a static CZM based on Turon et al. [33] is chosen in this study. This formulation was chosen for the reason described hereafter. For graphite/epoxy materials like the one used in this study, the fracture process zone (FPZ) is sufficiently small, indicating that LEFM is applicable in this case [33]. This implies that delamination growth, even under variable fracture mode loading conditions should be driven by fracture toughness and interlaminar strength should not affect the energy dissipated. But in a study, Turon et al. [34] shows that for existing CZM formulations, such as the one implemented in commercial codes like Abaqus, the interlaminar strengths can influence energy dissipation in variable mode loading scenarios such as the mixed mode bending. The formulation utilised in this study, addresses these inaccuracies in predicting energy dissipation under variable mode loading conditions. It has been demonstrated in Turon et al. [33] that this issue can be resolved by ensuring non-negative energy dissipation (i.e., preventing material healing) when the local fracture mode ratio changes during damage evolution. This introduced a new condition for the material properties of the cohesive model, as shown below.

$$ \frac{K_{sh}}{K_{33}} = \frac{G_{Ic}}{G_{IIc}} \left( \frac{\tau_{sh}^{0}}{\tau_{3}^{0}} \right)^{2} \quad (1) $$

where $K_{sh}$ and $K_{33}$ represent penalty stiffnesses, $\tau_{sh}^{0}$ and $\tau_{3}^{0}$ denote interlaminar strengths and $G_{Ic}$ and $G_{IIc}$ are critical energy release rates, corresponding to shearing and opening mode respectively. It should be noted that for graphite/epoxy materials like the one used in this study, mode II (shearing) and mode III (tearing) parameters are assumed to be the same. In the expression above, strengths and critical energy release rates are physical parameters determined experimentally, whereas penalty stiffnesses are numerical parameters. For the simulations the value of $K_{sh}$ is determined from the above expression. In general, the model incorporates mode-dependent penalty stiffness to ensure accurate simulation results under varying loading conditions. The computational algorithm of the damage model is implemented as a user-defined material subroutine (UMAT) in Abaqus/Implicit. The implementation of the UMAT is outlined in the flowchart in Fig. A.20.

The CZM damage model is implemented for cohesive elements of type COH3D8, which are three-dimensional elements with eight nodes and four integration points. Each element consists of two surfaces connected to the adjacent plies, initially separated by a near zero
<figure><img src="image_3.png" /><figcaption>Fig. 1. Details of the geometry, layup and delamination locations selected in this study.</figcaption></figure>

<table><tr><td>Layup</td><td>[0<sub>3</sub>/-45<sub>3</sub>/45<sub>3</sub>/90<sub>3</sub>]<sub>s</sub></td></tr><tr><td>Compressive force (F)</td><td>8 kN</td></tr><tr><td>Laminate length (L)</td><td>100 mm</td></tr><tr><td>Delamination length (l)</td><td>Starting with 10 mm</td></tr><tr><td>Laminate width (b)</td><td>20 mm</td></tr><tr><td>Laminate thickness (t<sub>Lam</sub>)</td><td>3.144 mm</td></tr><tr><td>Ply thickness (t<sub>Ply</sub>)</td><td>0.131 mm</td></tr></table>

distance. When subjected to loading, these surfaces separate from each other to facilitate the delamination growth. The initial thickness between two surfaces has to be small, since the constitutive response is modelled using a traction separation description. According to the Abaqus documentation [35], a thickness of 0.0001 mm or greater is recommended to avoid numerical problems, which led to the decision of using the value of 0.0001. The relationship between traction and separation of the two surfaces follows the implemented constitutive law. The traction separation law is given by

$$ \tau_i = (1 - d) [K_{ij}\Delta_j + \delta_{i3}K_{33}(-\Delta_3)] - \delta_{i3}K_{33}(-\Delta_3) \quad (i, j = 1..3) \quad (2) $$

where the constitutive tensor is uncoupled and diagonal with three stiffness terms $K_{33}$, $K_{11}$ and $K_{22}$ corresponding respectively to opening, shearing and tearing delamination modes, where the shearing and tearing mode stiffnesses are assumed to be same ($K_{sh}$) [33]. The variables $\Delta_j$ and $d$ represent the displacement jump corresponding to individual modes and the damage variable. It is important to note that the first term of the interface traction with the damage variable serves as the cohesive term for the opening mode, while the second term without damage acts as the contact term, valid under compression. The term with $\delta_{i3}$ represents the Kronecker delta, and $\langle \rangle$ denotes Macaulay brackets defined as $\langle x \rangle = \max(0, x)$.

As illustrated in Fig. A.20, the initial two inputs from ABAQUS to the user material (UMAT) are STRAN(3) and DSTRAN(3), representing the three modes of total displacement jump at the start of the increment and their increment during that specific time increment. Additionally, ABAQUS provides the time increment (DTIME), which is later used in the calculation of the viscous damage variable at the current time, $d'_v$. The damage variable and the viscous damage variable are passed on as state variables in UMAT, ensuring that their values from the previous time increment are available as inputs. The user also inputs the material parameters, which are listed in the flowchart (Fig. A.20).

The next step in Fig. A.20, shows the formulation based on [33]. An important parameter to define is the damage threshold $r'^{-1}$ to test whether the equivalent displacement jump $\Delta^o_{eq}$ in the current increment is large enough to cause damage. If $\Delta^o_{eq} > r'^{-1}$, then the new damage variable is calculated based on the linear evolution law shown below. Otherwise, the damage variable retains the same value as the previous increment.

$$ d^t = \frac{\Delta^f_{eq} (\Delta_{eq} - \Delta^o_{eq})}{\Delta_{eq} (\Delta^f_{eq} - \Delta^o_{eq})} \quad (3) $$

In the practical implementation of the CZM, it is common to encounter convergence difficulties, which can be overcome using a viscous regularisation technique based on Duvaut-Lions approach for rate-dependent viscoelastic materials [36]. In this technique, the damage variable $d'$ at the current time increment is replaced by the viscous

damage $d'_v$ in the current increment, whose time rate is proportional to $d' - d'_v$ with a numerical viscosity coefficient $\mu$ having units of time.

$$ \dot{d}_v = \frac{1}{\mu} (d' - d'_v) \quad (4) $$

In the calculation of tractions in Eq. (2), the damage variable $d'$ is replaced by $d'_v$ when the viscous regularisation technique is used. In the present study, a small enough viscosity coefficient has been applied to ensure it does not significantly affect the final solution. Additionally, besides the updated tractions, the UMAT must deliver a proper consistent tangent to the implicit solver every time the UMAT is called at an integration point. While the consistent tangent does not directly influence the simulation results, its quality is crucial for the convergence rate within the iterative Newton-Raphson scheme used by the implicit solver. Typically, an analytical expression is sought for the constitutive tangent [37]. However, such an expression can become overly complex for constitutive laws like the one used in this study. To address this challenge, the constitutive tangent tensor is computed numerically using a perturbation scheme [38].<fnref n="1" />

## 2.3. Structural modelling aspects

Fig. 2 shows the representative deformed configuration of the laminate, where the initial delamination and the resulting upper sublaminate buckling can be visualised. In Fig. 3, a schematic representation of the steps followed in the FE analysis is illustrated.

Initially, the details of the FE model have been illustrated, involving the selection of element type, mesh convergence, material parameters, and other relevant factors as shown in Fig. 3. The laminate has been modelled using incompatible mode eight-node brick elements (C3D8I) with incompatible deformation modes that eliminate the parasitic shear locking observed under bending for linear fully integrated solid elements [35]. This linear element has been chosen to ensure that the nodal locations are compatible with the adjacent linear eight-noded cohesive element COH3D8 for modelling delamination growth. The cohesive elements are incorporated along the horizontal plane passing through the initial delamination (see Fig. 2).

The mesh characteristics are selected based on mesh convergence studies of both postbuckling out-of-plane deformation and delamination growth behaviour which is based on the studies on cohesive zone length [30,39]. The mesh sensitivity studies has shown that modelling each three ply sublaminate (eg: $0_3$, see layup details in Fig. 1) with a single element through the thickness and in-plane dimension of 0.5 mm is sufficient to produce a compromise of accuracy and calculation

<footnote n="1">¹ For the detailed explanation of the background of CZM, the reader may refer to Turon et al. [33].</footnote>
<figure><img src="image_4.png" /><figcaption>Fig. 2. FE model showing details of the mesh and boundary conditions.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 3. Flowchart of the FE framework.</figcaption></figure>

For the accurate energy dissipation during delamination growth the mesh size along the loading direction must be sufficiently fine to ensure that enough interface elements exist within the fully developed cohesive zone length near the delamination tip. The length of this fully developed cohesive zone in a deformed geometry is the distance between the delamination tip where the complete damage has been reached and the point ahead of the tip at which the maximum cohesive traction is attained [39]. In order to find the appropriate mesh size a convergence study has been performed using the mode I double cantilever beam (DCB) test. This test has been selected because, for typical graphite/epoxy composite materials, the mode I cohesive zone length is the shortest and thus dictates the requirement for mesh size [30]. In Fig. 4 geometric details of the DCB model are presented

along with the comparison of the FE force-displacement with the analytical solution based on corrected beam theory and LEFM which relates mode I energy release rate to change of compliance with respect to change of crack length [40]. In the delamination growth region for a mesh size of 0.1 mm along the growth direction and 0.5 mm along the width direction, the DCB results converged to the analytical solution mentioned above. The mesh sizes that resulted from the DCB test have been used for modelling the fine mesh regions in the laminate for buckling analysis (see Fig. 2). However, only the plies adjacent to the cohesive elements have been modelled with fine mesh and the remaining plies have been modelled coarse with a size of 0.5 mm as discussed above.
<figure><img src="image_6.png" /><figcaption>Fig. 4. (a) FE model of DCB test and (b) Comparison of post mesh convergence FE solution with the analytical solution.</figcaption></figure>

<caption>**Table 1**</caption>

Linear-elastic and transversely isotropic ply properties of IM7/8552 unidirectional graphite/epoxy composite.

Source: Taken from Joosten et al. [41].

<table><caption>**Table 2**</caption><thead><tr><th colspan="6">Ply elastic properties</th></tr><tr><th>E<sub>11</sub> [MPa]</th><th>E<sub>22</sub> = E<sub>33</sub> [MPa]</th><th>v<sub>12</sub> = v<sub>13</sub> [-]</th><th>v<sub>23</sub> [-]</th><th>G<sub>12</sub> = G<sub>13</sub> [MPa]</th><th>G<sub>23</sub> [MPa]</th></tr></thead><tbody><tr><td>160 580</td><td>11 500</td><td>0.31</td><td>0.49</td><td>5430</td><td>3840</td></tr></tbody></table>

Interface cohesive zone properties of IM7/8552 unidirectional graphite/epoxy composite.

Source: Taken from Joosten et al. [41].

<table><thead><tr><th colspan="5">Cohesive zone properties</th></tr><tr><th>r<sub>1</sub><sup>o</sup> [MPa]</th><th>r<sub>2</sub><sup>o</sup> = r<sub>2h</sub><sup>o</sup> [MPa]</th><th>G<sub>1c</sub> [N/mm]</th><th>G<sub>11c</sub> = G<sub>11c</sub> [N/mm]</th><th>η<sub>BK</sub></th></tr></thead><tbody><tr><td>62.3</td><td>92.3</td><td>0.21</td><td>0.77</td><td>2.1</td></tr></tbody></table>

The linear-elastic and transversely isotropic properties of the unidirectional graphite/ epoxy composite ply (IM7/8552) have been used for the plies in this study and are listed in Table 1. Additionally, the cohesive zone properties of the lamina interface are given in Table 2.

The next step in modelling concerns defining the boundary conditions and interactions (see Fig. 3). To simulate compression, all displacement degrees of freedom at the left end of the laminate have been constrained and the translations on the other face are coupled to a reference point to which the compression force has been applied (see Fig. 2). TIE constraints are employed to constrain the active degrees of freedom on one surface to be same as the adjacent surfaces (ply-ply and ply-cohesive zone). To prevent penetration of adjacent surfaces at the location of the initial delamination, a surface-to-surface contact formulation has been used with default 'hard contact pressure-overclosure formulation'. The surface-to-surface discretisation offers advantage compared to traditional node-to-surface discretisation. In node-to-surface discretisation each slave node on one side of a contact interface effectively interacts with a point of projection on the master surface on opposite side of contact interface. Thus, the slave nodes

are constrained not to penetrate into the master surface, however, the nodes of the master surface can in principle penetrate into the slave surface. The surface-to-surface formulation solves this by enforcing contact conditions in an average sense over regions nearby slave nodes rather than only at the nodes. The pressure overclosure formulation defines how much contact pressure is created between two contact surfaces in the normal direction for a certain amount of penetration between them. The default hard contact formulation is physically most accurate which does not allow any penetration. 'Frictionless tangential behaviour' is assumed, leading to absence of tangential forces between the contact surfaces [35].

In the FE model of a delaminated laminate, buckling cannot be initiated in a flat configuration without any imperfection [2]. For this reason, similar studies in the literature [2,5] have introduced geometric imperfection. Following this approach, this study incorporates geometric imperfection based on the LBA as illustrated in the second step of Fig. 3. The first eigenmode extracted from the LBA has been chosen to define the imperfection because of its plausible shape and the lowest magnitude of eigenvalue. For the nonlinear analysis, the FE model is seeded with an imperfection calculated by multiplying the deformation of the first eigenmode by a user-specified factor. For instance, if the factor is 0.01, this means that 1% of the calculated normalised deformation (maximum value of 1.0) of the first eigenmode is applied to the undeformed FE mesh. In the studies examining post-buckling behaviour of delaminated composite laminates in Mekonnen et al. [2] and Mohammadi and Shahabi [5], imperfection magnitudes ranging from 0.4% to 5% have been generally employed. In another study by Raimondo et al. [42] investigating delamination growth in postbuckled stiffened composite panel an imperfection magnitude of 1% of the first eigenmode has been utilised. In the present study, the smaller imperfection magnitudes of 0.1% and 1.0% have been compared to demonstrate their influence on the results and thereby leading to important conclusions discussed later.

Final step involved performing a geometric nonlinear analysis (GNA) using Full Newton solution technique (see Fig. 3). Since deformations are large in buckling analysis, geometric nonlinearity has been accounted for by enabling the NLGEOM in ABAQUS [35].

## 2.4. Validation of the FE model

To validate the modelling aspects presented in Section 2.3 and the CZM detailed in Section 2.2, the numerical results documented by Wang and Zhang [18] and the corresponding experimental results from Kutlu and Chang [26] have been utilised. These results have also served as benchmark for other studies in literature on delamination buckling [2,5,19,20]. In this example, the buckling response of a laminate featuring a single through-width delamination has been studied. The details of geometry and layup of the laminate has been presented in Fig. 5(a) while the material is Fiberite T300/976 graphite/epoxy, with its parameters listed in Tables 3 and 4. In the experiment in Kutlu and Chang [26], the artificial delamination was created by introducing a Teflon film of thickness 0.0254 mm. This delamination created by Teflon film divides the laminate into thin upper and thick lower sublaminates. The effect of Teflon film can be accounted for in simulation by creating a corresponding gap and causing an undulation of thin sublamine. However, in this study, the effect of the Teflon film has not been explicitly modelled. Instead, the delaminated laminate is modelled without a gap, and the geometric imperfection is introduced using LBA to initiate buckling. This approach is consistent with the numerical model of Wang and Zhang [18], which served as a reference. Similarly, other studies in the literature [2,5] also do not explicitly account for the gap introduced by the Teflon film. Furthermore, the gap created by the Teflon film may lead to resin accumulation at the delamination front, potentially increasing resistance to growth. However, there is no mention of this effect in the selected literature [18,26], suggesting that its influence may have been minimal.
<caption>**Table 3**</caption>

Linear-elastic and transversely isotropic ply properties of Fiberite T300/976 graphite/epoxy composite.

Source: Taken from Wang and Zhang [18].

<table><caption>**Table 4**</caption><thead><tr><th colspan="7">Ply elastic properties</th></tr><tr><th>E<sub>11</sub> [MPa]</th><th>E<sub>22</sub> = E<sub>33</sub> [MPa]</th><th>ν<sub>12</sub> = ν<sub>13</sub> [-]</th><th>ν<sub>23</sub> [-]</th><th>G<sub>12</sub> = G<sub>13</sub> [MPa]</th><th>G<sub>23</sub> [MPa]</th></tr></thead><tbody><tr><td>139 300</td><td>9720</td><td>0.29</td><td>0.40</td><td>5580</td><td>3450</td></tr></tbody></table>

Interface cohesive zone properties of Fiberite T300/976 graphite/epoxy composite.

Source: Taken from Wang and Zhang [18] and Mekonnen et al. [2]

<table><thead><tr><th colspan="5">Cohesive zone properties</th></tr><tr><th>τ<sup>a</sup><sub>1</sub> [MPa]</th><th>τ<sup>a</sup><sub>2</sub> = τ<sup>a</sup><sub>3</sub> = τ<sup>a</sup><sub>4</sub> [MPa]</th><th>G<sub>1c</sub> [N/mm]</th><th>G<sub>11c</sub> = G<sub>111c</sub> [N/mm]</th><th>η<sub>BK</sub></th></tr></thead><tbody><tr><td>44.54</td><td>106.90</td><td>0.0876</td><td>0.3152</td><td>1.45</td></tr></tbody></table>

<figure><img src="image_7.png" /><figcaption>(a)</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>(b)</figcaption></figure>

**Fig. 5.** (a) Comparison of compressive force vs. displacement curve between Wang and Zhang [18] and present study (b) Comparison of load deflection history between Wang and Zhang [18] and present study (The schematic representation of the investigated laminate and its geometric parameters is included in the figure, where // designates the through-thickness location of delamination).

The delaminated laminate has been subjected to gradually increasing compressive loading until the laminate failure is observed. Comparison of the results with the numerical investigation based on the B-spline finite strip method of Wang and Zhang [18] is illustrated in Fig. 5.

Fig. 5(a) depicts the axial compression force versus displacement curve while Fig. 5(b) shows the graph of axial compression force versus out-of-plane displacement at the centre of the upper and lower laminate surfaces. As mentioned above Fig. 5(a) also depicts details of the laminate, including the layup and its dimensions. A further comparison has been made for the same problem with the numerical results reported by Wang and Zhang [18] and the experimental results in Kutlu and Chang [26], as illustrated in Fig. 6. Fig. 6(a) and (b) showcase the axial compression force versus local axial strains at the centre of the upper and lower laminate surfaces, with the specific locations highlighted in the figures.

<figure><img src="image_9.png" /><figcaption>(a)</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>(b)</figcaption></figure>

**Fig. 6.** Comparison of axial strain histories of the present study with the experimental results of Kutlu and Chang [26] and the numerical predictions of Wang and Zhang [18] on the middle of (a) upper sublamine surface and (b) lower sublamine surface.

The results of the current model have been plotted for imperfection magnitudes of 0.1% and 1.0% and it reveals that the main difference lies in the buckling behaviour of thin upper sublamine (Fig. 5(b)). It is observed that the upper sublamine deviates from the primary compression path at lower loads compared to the findings of Wang and Zhang [18]. The results indicate that the initial upper sublamine buckling occurs approximately at a load of 0.93 kN for 0.1% and 0.80 kN for 1.0% in the present study, aligning closely with literature [18], which reported a value of 1.015 kN. The slight deviation is attributed to a geometric imperfection introduced in the present model for nonlinear analysis.

Examining the experimental data of local surface strains in Fig. 6(a) [26], a snap-through from compression to tension strain occurs at a higher load of 1.39 kN compared to both the present study and the numerical study by Wang and Zhang [18]. This transition is caused by the buckling of the upper sublamine. In Kutlu and Chang [26], the delayed occurrence of upper sublamine buckling in the experimental results has been attributed to the sticking of the upper sublamine to the lower sublamine through the Teflon film. This buckling of the upper sublamine results in a reduction of axial stiffness, evident in the change of slope beyond 1 kN in the compressive force versus displacement curve in Fig. 5(a). At load levels of 1.33 kN in [18] and 1.30 kN in the present study (for both imperfection magnitudes), the laminate undergoes an initial unstable delamination growth. In both the present study and literature [18], the delamination extends initially by around 10 mm, marked by a distinctive load drop. Subsequently, with increasing load, the delamination gradually grows towards the two ends of the laminate. After the delamination has grown to the two ends, the final global buckling occurs. This global buckling corresponds to the buckling of the thicker lower sublamine on the opposite side at a load level of around 2.32 kN for both the cases, as observed in Fig. 5(b). The buckling of thicker sublaminates at the higher load is explained by its superior bending stiffness. The surface strain histories presented in Fig. 6 overall compare well with the experimental results of Kutlu
<figure><img src="image_11.png" /><figcaption>Fig. 7. For near surface delamination at imperfection magnitude of 0.1%: (a) Contour plot of evolution of out-of-plane displacement and (b) Damage development in cohesive elements with normalised delamination length. For enhanced visualisation the magnification scale of out-of-plane displacement plots have been increased to 3.0.</figcaption></figure>

and Chang [26] and the numerical results of Wang and Zhang [18]. The global buckling failure in the experiment of Kutlu and Chang [26] occurred however at a slightly higher load of 2.44 KN. In summary, the results from the proposed modelling strategy are in good agreement with existing studies in the literature and can therefore be used for further investigation. Wang and Zhang [18] and several other studies [2,5] in the literature have primarily investigated the force-displacement response of laminates with through-width delaminations. However, these studies do not explore the factors influencing the unstable growth, which are decisive for implementing the 'slow growth' design concept and are therefore the main objectives of the current study.

## 3. Results and discussions

Following the modelling strategy detailed in Sections 2.2 and 2.3, the parametric study outlined in Section 2.1 has been conducted and this section provides a comprehensive discussion of the obtained results. The initial discussion focuses on the influence of through-thickness delamination location in the selected QI laminate. This study investigates the occurrence of unstable growth in the QI laminate for each through-thickness location of the delamination, identifying the most critical configuration. For the most critical case in this QI laminate, the influence of laminate dimensions has been investigated. Subsequently, the original laminate is compared to a stiffer one to gain deeper insights into the problem.

### 3.1. Influence of through-thickness delamination location in a QI laminate

#### 3.1.1. Near surface delamination

In this section the results of a near surface delamination defect are presented. Figs. 7 and 8 are plotted at the reference load of 8 KN. In Fig. 7, the left side displays contour plots of the out-of-plane displacement, w, while the right side illustrates damage growth in cohesive elements

using damage variable, d across various delamination lengths. The plots in Fig. 7 have been presented upto fine mesh region of the laminate, which extends 15 mm on either side of the delamination. Fig. 8(a) shows the graph of the out-of-plane displacement, w, at the centre of the upper and lower laminate surfaces in relation to delamination lengths, while Fig. 8(b) illustrates the change in the axial stiffness of the laminate with respect to delamination length. Axial stiffness here refers to the tangent stiffness derived from the force-displacement relationship which, has been elaborated in the Appendix B. In all the figures the delamination length, *l* has been normalised with respect to the laminate length *L*. Initially, simulations have been conducted by varying the normalised delamination length in increments of 0.1. However, to gain a deeper understanding of the problem, further analyses have been performed with finer increments, as shown in the insets of Fig. 8.

Fig. 7 just reports the contour plot for imperfection magnitude of 0.1% whereas in Fig. 8 a comparison of different imperfection magnitudes are presented. The initial delamination divides the laminate into two sublaminates, with the upper slender sublamine possessing significantly lower bending stiffness compared to the lower thick sub-lamine. Initially the results corresponding to imperfection magnitude of 0.1% have been presented. For this imperfection magnitude, no significant out-of-plane displacement of both the sublaminates has been observed until the normalised delamination length of 0.13 (see Fig. 8(a)). In the global graph of Fig. 8(a), when the normalised delamination length is increased from 0.13 to 0.14, the out-of-plane displacement of the upper sublamine increases to 0.0988 mm, exceeding 25% of its thickness (0.393 mm). This significant increase indicates that the upper sublamine has undergone buckling and the laminate is be said to have undergone local buckling. However, this global graph does not provide a detailed visualisation of the buckling behaviour. Therefore, additional simulations have been conducted for intermediate normalised delamination lengths of 0.1325, 0.1350, and
<figure><img src="image_12.png" /><figcaption>Fig. 8. For near surface delaminations: (a) Evolution of out-of-plane displacement with normalised delamination length and (b) Axial stiffness change with normalised delamination length.</figcaption></figure>

0.1375. The results of these have been shown in the inset in Fig. 8(a). It is shown that unlike the assumed linear increase in out-of-plane displacement shown in the global graph, the actual behaviour is gradual and nonlinear. This indicates that local buckling does not occur at a discrete point but rather as a continuous process, with out-of-plane displacement gradually increasing as the delamination length increases.

As shown in Fig. 7(b), for normalised delamination length of 0.14, local buckling causes some damage at the edges of the delamination front. This combined effect of local buckling and damage results in a 4.5% reduction in axial stiffness compared to the original value (see Fig. 8(b)). In addition to global graph, again an inset has been shown in Fig. 8(b) which illustrates the decrease in stiffness for intermediate normalised delaminations between 0.13 and 0.14. The inset shows a more gradual stiffness decrease compared to the global graph, where a linear reduction is assumed. This behaviour is expected and aligns with the out-of-plane displacement growth shown in the inset of Fig. 8(a). With increase in normalised delamination length to 0.15 the out-of-plane displacement increases further beyond the point of local buckling. For upper sublamine it increases to 0.254 mm and for lower sublamine it increases marginally to 0.053 mm. This causes damage to develop along the entire delamination front (see Fig. 7(b)) and leads to a 15% reduction in axial stiffness compared to the original value (see Fig. 8(b)). At normalised delamination length of 0.16, the delamination extends in an unstable manner by 22 mm (see Fig. 7(b)), causing out-of-plane displacement in the lower sublamine also to increase significantly (see Fig. 8(a)). At this point, the laminate is considered to have undergone mixed local/global buckling failure. In Fig. 8(b) this is indicated by sudden drop in the value of axial stiffness.

Introducing a higher imperfection magnitude of 1.0% in the lamina results in the thin sublamine exhibiting increased out-of-plane displacement even at smaller delamination lengths. (see Fig. 8(a)).

In this case, the out-of-plane displacement gradually begins to increase after the normalised delamination length exceeds 0.11. It remains higher than the values for the imperfection magnitude of 0.1% until the normalised delamination length reaches 0.14. However, as the normalised delamination length is increased to 0.15, the out-of-plane displacement values nearly match between the two imperfection magnitudes. Likewise, at a normalised delamination length of 0.16, the delamination extends in an unstable manner by 22 mm, resulting in laminate failure due to mixed local/global buckling. This outcome is similar to the trend observed in the validation example in Section 2.4, where although the initial buckling behaviour of the upper thin sublamine differed, the delamination growth load did not differ much for the considered imperfection magnitudes. Therefore, the results for both imperfection magnitudes show that unstable growth occurs at a critical delamination length of 16 mm, and the delamination extends by 22 mm. The question arises as to what factors influence this critical delamination size and whether unstable growth can be prevented, both of which are addressed later.

### 3.1.2. Intermediate delamination

This section presents the results of intermediate delamination defect. Fig. 9(a) presents the graph of the evolution of out-of-plane displacement at the centre of laminate surfaces with respect to the normalised delamination length, while Fig. 9(b) shows the change in axial tangent stiffness in relation to the normalised delamination length. In Fig. 10 FE contour plots have been presented corresponding to the critical delamination length at which the unstable growth and laminate failure occurs. Fig. 10(a) illustrates the out-of-plane displacement while Fig. 10(b) depicts the delamination growth with the help of damage variable, d. Unlike the near-surface delamination case, the contour plots
<figure><img src="image_13.png" /><figcaption>Fig. 9. For intermediate delaminations: (a) Evolution of out-of-plane displacement with normalised delamination length and (b) Axial stiffness change with normalised delamination length.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Fig. 10. For intermediate delamination: (a) Out-of-plane displacement and (b) Damage growth by more than 30 mm in cohesive elements at the normalised delamination length of 0.22. Magnification scale of 3.0 is used for out-of-plane displacement plot to keep it consistent with other plots.</figcaption></figure>

are shown for a region just beyond the fine mesh area of the laminate, which extends 15 mm on either side of the delamination.

Similar to near surface configuration, the initial delamination divides the laminate into two sublaminates ([0<sub>3</sub>/−45<sub>3</sub>] and [45<sub>3</sub>/90<sub>3</sub>/90<sub>3</sub>/45<sub>3</sub>/−45<sub>3</sub>/0<sub>3</sub>]), with the upper slender sublaminate possessing lower bending stiffness compared to lower thick sublaminate. For this reason, the progression of buckling follows the same trend as near-surface delamination, beginning with local buckling and transitioning to mixed local/global buckling failure as the delamination length is increased. For an imperfection magnitude of 0.1%, the out-of-plane displacement of the thin sublaminate gradually increases to 0.0714 mm at a normalised delamination length of 0.19, which corresponds to 9.0% of the thickness (which is 0.786 mm) of the thin sublaminate (see Fig. 9(a)). As the normalised delamination length increases from 0.19 to 0.20, the out-of-plane displacement of the thin sublaminate significantly increases from 9.0% to 15.0% of its thickness, indicating that the laminate has undergone local buckling. Up to a normalised delamination length of 0.19, the gradual increase in the out-of-plane displacement of the thin sublaminate can be attributed to bending due to assymetry.

Compared to near surface delamination defect, local buckling occurs for higher delamination length due to the greater bending stiffness of the thin sublaminate. At a normalised delamination length of 0.20,

local buckling leads to a 3.8% reduction in axial stiffness compared to the axial stiffness at a normalised delamination length of 0.10 (see Fig. 9(b)). As the normalised delamination length increases to 0.21, the out-of-plane displacement continues to increase beyond the point of local buckling. For the upper sublaminate, it increases to 0.236 mm, while for the lower sublaminate, it marginally increases to 0.0488 mm (see Fig. 9(a)). At a normalised delamination length of 0.22, the delamination extends in an unstable manner by more than 30 mm which is beyond the fine mesh region. However, the results indicate that even as the delamination grows to the edges of the fine mesh region, there is a notable increase in the out-of-plane displacement in the lower sublaminate (see Fig. 10(a)). Similar to a near-surface delamination configuration, the laminate at this point has experienced mixed local/global buckling failure. Consequently there is sudden drop in the axial stiffness as shown in Fig. 9(b).

Increasing the imperfection magnitude to 1.0% does not notably influence the results. The pre-local buckling behaviour shows slightly higher out-of-plane displacements in the thin sublaminate compared to a lower imperfection magnitude (see Fig. 9(a)). This results in further reduction of axial stiffness as compared to imperfection magnitude of 0.1% (see Fig. 9(b)). However, the normalised delamination length, corresponding to unstable delamination growth, remains unchanged at
<figure><img src="image_15.png" /><figcaption>Fig. 11. For midplane delamination: Evolution of out-of-plane displacement with normalised delamination length. Magnification scale of 3.0 is used to keep it consistent with other plots.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Fig. 12. For midplane delamination: Axial stiffness change with normalised delamination length.</figcaption></figure>

0.22. This outcome is similar to the near surface configuration. Similar to the near-surface configuration, unstable growth is also observed in this intermediate configuration, but at a higher critical delamination length of 22 mm. As the qualitative behaviour in this case is similar to the near-surface configuration, only the near surface configuration has been investigated further.

### 3.1.3. Midplane delamination

For midplane delamination, in Fig. 11, evolution of out-of-plane displacement has been shown in relation to the normalised delamination length. It should be noted that these contour plots display entire length of the laminate. Fig. 12 illustrates the change of axial tangent stiffness with the normalised delamination length. Fig. 11 displays the displacements corresponding to an imperfection magnitude of 1.0%, while Fig. 12 compares the axial stiffness change for both imperfection magnitudes used in this study. The midplane delamination divides the laminate into two sublaminates with identical bending stiffness. Consequently the laminate does not undergo local buckling initially as demonstrated for previous configurations. Instead both the

sublaminates tend to deform in the same direction as shown in Fig. 11. The deformation increases gradually as the delamination length is increased upto a normalised delamination length of 0.40. At a normalised delamination length of 0.45, the laminate exhibits significant out-of-plane displacement with U-shaped deformation. This indicates that the laminate has undergone failure in global buckling mode.

Fig. 12 shows that increasing the normalised delamination length from 0.10 to 0.25 results in only a marginal reduction in axial stiffness, with decreases of 2.9% for an imperfection magnitude of 0.1% and 4.3% for an imperfection magnitude of 1.0%. This differs from the near-surface and intermediate delamination configurations, where the laminate failed at the lower normalised delamination lengths of 0.16 and 0.22, respectively. This shows that mode of buckling has significant influence on the behaviour of the laminate. Beyond the normalised delamination length of 0.25, the axial stiffness continues to decrease with the delamination length, until reaching a normalised length of 0.40. The subsequent simulation at a normalised delamination length of 0.45 exhibits a sudden drop in axial stiffness, indicating global buckling failure as previously discussed.

The simulation results also show that, in contrast to previous configurations, no unstable growth of delamination is observed in this configuration for the same reference load of 8 kN. However, if the applied load is increased marginally to 9.2 kN, then unstable growth and failure of laminate is observed at the critical delamination size of 30 mm, which suggests that for this midplane delamination, the load level can be optimised to prevent unstable growth. The next sections similarly investigates the most critical near-surface configuration, to identify various factors influencing the unstable growth behaviour like the laminate dimensions and sublamine stiffness ratios.

### 3.2. Influence of laminate dimensions

Initially, for the near-surface delamination configuration, laminates of two lengths, *L*, measuring 50 mm and 100 mm, are compared. The Fig. 13 below presents a comparison between the two laminate lengths. Fig. 13(a) illustrates the evolution of out-of-plane displacements for the
<figure><img src="image_17.png" /><figcaption>**Fig. 13.** Comparison between different laminate lengths for near surface delamination (a) Evolution of out-of-plane displacement with delamination length and (b) Axial stiffness change with delamination length.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>**Fig. 14.** For near surface delamination in laminate with length of 50 mm (a) Out-of-plane displacement and (b) Damage growth by 6 mm in cohesive elements at delamination length of 18 mm. Magnification scale of 3.0 is used for out-of-plane displacement to keep it consistent with other plots.</figcaption></figure>

thin and thick sublaminates as a function of delamination length, while Fig. 13(b) shows the reduction in axial tangent stiffness with increasing delamination length.

Note that the delamination length in the figures has not been normalised with respect to the laminate length, as the laminate lengths differ in the two cases discussed, making direct comparison not possible. The original length, *L* = 100 mm, has already been discussed in Section 3.1.1. In this case, once the delamination length exceeds 13 mm, the out-of-plane displacement in the thin sublaminate begins to increase, indicating the onset of local buckling. Beyond this point, the out-of-plane displacement continues to increase with the increasing delamination length. At a delamination length of 15 mm, the out-of-plane displacement in the thin sublaminate reaches 0.254 mm, while the thick sublaminate experiences a marginal increase to 0.053 mm. When the delamination length is increased to 16 mm, there is unstable growth of delamination by 22 mm and a significant increase in out-of-plane displacements in both the thin and thick sublaminates. For the short laminate length of 50 mm, the out-of-plane displacement in the thin sublaminate remains identical to that observed for laminate length of 100 mm up to a delamination length of 15 mm. This occurs

because, under the same applied load, the buckling behaviour of the thin sublaminate is governed by the delamination length, which defines the free length. Since the free length is identical in both cases, the behaviour remains the same. However, the out-of-plane displacement of the thick sublaminate differs between the two cases, as it is influenced by the overall laminate length. For the original laminate length of 100 mm, the higher out-of-plane displacements of thick sublaminate leads to unstable delamination growth by 22 mm at a delamination length of 16 mm. In contrast, for the short laminate length, the unstable growth occurs at a delamination length of 18 mm, with the delamination extending by only 6 mm in this case (see Fig. 14). It is also observed that, even at the point of unstable growth, the out-of-plane displacement of the thick sublaminate does not increase significantly for the short laminate length. However, the impact of unstable growth is evident in Fig. 13(b), where a sudden decrease in axial tangent stiffness is observed for 18 mm delamination length. Therefore, it is demonstrated that varying the laminate length can affect the growth characteristics, primarily due to changes in the buckling behaviour of the thick sublaminate.
<figure><img src="image_19.png" /><figcaption>**Fig. 15.** Comparison between different laminate widths *b* for near surface delamination (a) Evolution of out-of-plane displacement with normalised delamination length and (b) Axial stiffness change with normalised delamination length.</figcaption></figure>

Next, the effect of width is analysed by comparing a laminate with the original width, *b* of 20 mm to one with a width of 10 mm. Accordingly, the reference load is proportionally reduced from 8 kN to 4 kN for laminate with low width. Similar to the previous case, Fig. 15 compares for the two widths: (a) the evolution of out-of-plane displacements with normalised delamination length and (b) the reduction in axial tangent stiffness with normalised delamination length.

For the laminate with a 10 mm width, a trend similar to the original laminate is observed. As the normalised delamination length increases from 0.13 to 0.14, the out-of-plane displacement in the thin sublaminate increases significantly, indicating the onset of local buckling. However, at a normalised delamination length of 0.14, the out-of-plane displacement in the thin sublaminate is higher for the 10 mm width laminate, reaching 0.1446 mm, compared to 0.0988 mm for the original width laminate. Similarly, at a normalised delamination length of 0.15, the out-of-plane displacements (for thin and thick sublaminates) in the 10 mm wide laminate are higher than that of the original laminate. When the normalised delamination length is slightly increased to 0.1525, unstable delamination growth is observed, causing the delamination to extend by 25 mm (see Fig. 16). In the study of the original laminate until now, the normalised delamination length has been incremented in steps of 0.1, with unstable growth observed at a normalised delamination length of 0.16. However, it is possible that unstable growth could occur at a normalised length between 0.15 and 0.16. To investigate this, additional calculations have been performed at normalised delamination lengths of 0.1525, 0.1550, and 0.1575. Already at a normalised delamination length of 0.1575, unstable growth is observed in this case, causing the delamination to extend by approximately 22 mm. Thus, it is demonstrated that reducing the laminate width by half and proportionally reducing the reference

load can cause unstable growth to occur marginally earlier, while also resulting in a marginally greater delamination extension compared to the original laminate. However, its influence is not as significant as that of laminate length.

### 3.3. Comparison with stiffer laminates

The Section 3.1 has shown that the near-surface configuration is the most critical, with unstable growth and laminate failure occurring at a critical delamination length of 16 mm in a QI laminate with the layup [0<sub>3</sub>/−45<sub>3</sub>/45<sub>3</sub>/90<sub>3</sub>/90<sub>3</sub>/45<sub>2</sub>/−45<sub>3</sub>//0<sub>3</sub>], where ‘//’ indicates the location of the through-thickness delamination. However, in practice there can be laminates where the ratio of stiffness of thick to thin sublaminate can be higher. Therefore, in this section the influence of laminate stiffness is evaluated by considering a stiffer thick sublaminate. This is achieved by increasing the proportion of 0° plies in the thick sublaminate while maintaining the overall thickness and symmetry of the laminate. The considered 0-dominant laminate thus has a layup of [0<sub>3</sub>/−45<sub>2</sub>/0<sub>3</sub>/45<sub>2</sub>/90<sub>2</sub>/90<sub>2</sub>/45<sub>2</sub>/0<sub>3</sub>/−45<sub>2</sub>//0<sub>3</sub>] and the behaviour of two laminates have been compared at reference load of 8 kN (see Fig. 17).

Fig. 17(a) presents out-of-plane displacement with respect to normalised delamination length and Fig. 17(b) illustrates the change in axial tangent stiffness also in relation to normalised delamination length. For QI laminate it is already known that unstable growth occurs at critical delamination length of 16 mm and the delamination extends by 22 mm, resulting in laminate failure. On the other hand, for 0-dominant laminate the out-of-plane displacement of thin sublaminate shows notable increase (or local buckling) only beyond the normalised delamination length of 0.16. Beyond the normalised delamination length of 0.16 the out-of-plane displacement also begins to increase in thick sublaminate and this reflects in reduction of axial stiffness of the laminate. With further increase in the normalised delamination length, the out-of-plane displacement continues to increase and with that the axial stiffness also continues to reduce. However, for upto normalised delamination length of 0.40 no unstable growth of delamination is observed. This shows that onset of local buckling can be delayed and also the unstable growth can be prevented by increasing the stiffness of thick sublaminate. However, since local buckling does occur in the laminate beyond a normalised delamination length of 0.16, it can induce peel stresses and potentially lead to growth under fatigue conditions. Despite this, the growth is not expected to become unstable.

For the 0-dominant laminate, further investigation has been conducted to determine the load levels at which unstable growth occurs. The Fig. 18 below shows failure loads corresponding to various normalised delamination lengths in the zero-dominant laminate. For a normalised delamination length of 0.13, if the load is increased to its failure value of 15.9 kN, the delamination growth will be unstable, as the failure loads for subsequent delamination lengths are lower than 15.9 kN. A similar trend will be observed for load levels corresponding to the failure loads of normalised delamination lengths below 0.2. However, for normalised delamination lengths beyond 0.2, the load must be increased for further growth, as the failure load increases with increasing delamination length. The turning point in the variation of failure load after a normalised delamination length of 0.2 can be explained based on the energy release rate, *G*. A study by Chai et al. [43] employed an analytical model based on beam theory and a delamination growth criterion governed by the fracture mechanics-based energy release rate, *G*, offering some insights into the observed phenomena. It was shown that when a constant load is applied to the laminate while increasing the delamination length, *G* initially increases with the delamination length until a critical delamination length is reached. Beyond this point, *G* begins to decrease as the delamination length increases. It was shown that this reversal in the trend of *G* is responsible for the corresponding reversal in failure load beyond a critical value, which, in the present study, occurs at 0.2. However, to gain deeper insights into this phenomenon, more detailed analysis is
<figure><img src="image_20.png" /><figcaption>Fig. 16. For near surface delamination in laminate with width of 10 mm (a) Out-of-plane displacement and (b) Damage growth by 25.3 mm in cohesive elements at normalised delamination length of 0.1525. Magnification scale of 3.0 is used for out-of-plane displacement to keep it consistent with other plots.</figcaption></figure>

<figure><img src="image_21.png" /></figure>

<figure><img src="image_22.png" /><figcaption>Fig. 17. Comparison of QI laminate with 0-dominant laminate for near surface delamination (a) Evolution of out-of-plane displacement with normalised delamination length (b) Normalised axial stiffness change with normalised delamination length.</figcaption></figure>

needed using efficient and accurate computational models. Although, as explained earlier, this phenomenon has been reported earlier by Chai et al. [43], the problem was simplified by assuming the laminate to be isotropic and nonlinear stretching effects in the postbuckling regime was not considered. Alternatively, recent models, such as the one proposed by Köllner et al. [44], account for composite anisotropy and stacking sequence while also incorporating a nonlinear kinematic approach for the postbuckling regime. While these models can be used for gaining deeper insights, they still require further refinement. However, a detailed investigation of this aspect is not included within the scope of the current study. In general, at load levels exceeding 12 kN, the zero-dominant laminate is prone to unstable growth and the critical delamination length associated with unstable growth can be increased by decreasing the load. This demonstrates that for near-surface case, load levels can be optimised to delay or prevent unstable growth, which is a necessity for implementing the ‘slow growth’ concept. In the ‘no growth’ approach on the other hand, the stiffness of the thick

<figure><img src="image_23.png" /><figcaption>Fig. 18. Failure load corresponding to various normalised delamination length in 0-dominant laminate.</figcaption></figure>

sublaminate must be increased to further delay or prevent the onset of local buckling and the resulting peel stresses at the interface. This is discussed further for a laminate layup of [[0$_{3}$]/-45$_{3}$/45$_{3}$/90$_{3}$]$_{1}$$_{/}$$_{3}$/[90$_{3}$/45$_{3}$/-45$_{3}$/0$_{3}$]$_{1}$$_{/}$$_{2}$/90$_{3}$/45$_{3}$/-45$_{3}$/0$_{3}$], where thick sublaminate has significantly higher thickness and thus stiffness. The Fig. 19 below presents the graph of the evolution of out-of-plane displacement at the centre of laminate surfaces with respect to the normalised delamination length. The graph corresponds to same reference load of 8 kN. It is seen that compared to QI and zero-dominant laminate, the local buckling occurs much later at the normalised delamination length of 0.27. The thick sublaminate on the other hand shows negligible out-of-plane displacement due to much higher stiffness. Thus, for achieving no-growth condition number of plies may have to be increased in the thick sublaminate, leading to a weight penalty. To further investigate the growth characteristics, a parametric study involving various sublaminate stiffness ratios and load levels is necessary. Theoretical formulations, such as those by Köllner et al. [44], provide a foundation for conducting such investigations. However, this falls beyond the scope of the current study and can be considered as future research.

## 4. Conclusions and outlook

Delaminations in composite structures can be significantly detrimental to its mechanical performance by increasing susceptibility to premature buckling and further growth. This study explores whether the delamination growth can become unstable in the context of delamination buckling and identifies the key factors that govern the growth characteristics. Understanding the growth characteristics will aid designers in planning inspection schedules within the framework of damage tolerance concepts, such as the ‘slow-growth’ approach. This offers an advantage over applying the conservative ‘no growth’ concept, which can lead to weight penalties. A parametric ‘Effect of Defect’ study has been conducted, in which the applied load is held constant while the delamination length in the models is progressively increased. The objective is to evaluate whether the unstable growth could occur, and if so, determine the critical delamination length at which it starts. Initially
<figure><img src="image_24.png" /><figcaption>**Fig. 19.** Evolution of out-of-plane displacement with normalised delamination length for laminate configuration: [[0<sub>3</sub>/ - 45<sub>3</sub>/45<sub>3</sub>/90<sub>3</sub>]<sub>1</sub>/[90<sub>3</sub>/45<sub>3</sub>/ - 45<sub>3</sub>/0<sub>3</sub>]<sub>2</sub>/90<sub>3</sub>/45<sub>3</sub>/ - 45<sub>3</sub>//0<sub>3</sub>], where // designates the through thickness delamination location.</figcaption></figure>

the influence of through thickness location of delamination has been studied in a generic QI laminate. For delaminations near the surface and at intermediate through-thickness locations, the laminate initially experienced local buckling (thin sublaminate buckling), followed by unstable growth in both cases. The near-surface configuration is more critical of the two, due to earlier onset of local buckling. In the mid-plane delamination, the laminate shows a tendency for global buckling, and under the same applied load, no unstable growth is observed. This configuration is the least critical, with failure resulting from global buckling of the sublaminates at a significantly greater delamination length. However, by incrementally increasing the applied load, it is observed that the laminate can exhibit unstable growth even in the midplane configuration. This suggests that load levels can be optimised to prevent unstable growth. This study also evaluates the influence of laminate dimensions for most critical near surface configuration. The findings indicate that when the laminate length is reduced by half, the thick sublaminate exhibits reduced out-of-plane displacement, which in turn affects the critical delamination length at which unstable growth occurs and also the extent of growth. On the other hand, reducing the laminate width by half only has a marginal influence on the results. Finally a comparison has then been made between the previously discussed QI laminate and a stiffer, zero-dominant laminate for the near-surface configuration. At the same applied load, it is shown that no unstable growth occurs in the zero-dominant laminate. However, when the load is increased, unstable growth eventually emerges in the laminate, and with further increase in load, the critical delamination size decreases. It has also been shown that by significantly increasing the stiffness of thick sublaminate the local buckling and thus peel stresses at the interface can be significantly delayed. This would ensure that ‘no growth’ occurs at the smaller delamination lengths but it comes at the price of added weight penalty. To gain further insights into this problem, a more comprehensive investigation incorporating various sublaminate stiffness ratios and load levels can be carried out as the next step. This study lays the groundwork for identifying stable growth configurations in future research on the impact of delamination defects under fatigue loading. Given its critical nature resulting from local buckling, the near-surface configuration is considered the most relevant for fatigue studies.

CRediT authorship contribution statement

S.S. Venkat: Writing – original draft, Writing – review & editing, Conceptualization, Methodology, Investigation, Formal analysis, Visualization, Validation. S. Scheffler: Writing – review & editing, Methodology, Supervision, Resources, Project administration. P.M. Anilkumar: Writing – review & editing, Methodology, Visualization, Supervision, Project administration. E. Baranger: Writing – review & editing, Supervision, Resources, Funding acquisition. R. Rolfes: Writing – review & editing, Supervision, Resources, Funding acquisition.

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgements

The authors are grateful to the International Research Training Group 2657 (IRTG 2657) funded by the German Research Foundation (DFG) (Grant reference number 433082294). The authors also acknowledge the compute cluster of the Leibniz University Hannover, which has been mainly used to perform the presented simulations. The compute cluster is partly funded by DFG – Project-IDs INST 187/742-1 FUGG and INST 187/592-1 FUGG.

Appendix A. Implementation of static CZM

See Fig. A.20.

Appendix B. Axial stiffness of the laminate

This section describes the determination of the axial stiffness of the laminate, for which the force–displacement curve along with corresponding stiffness drop has been illustrated for the near-surface configuration with a delamination size of 16 mm (Fig. B.21). The 16 mm delamination represents the critical length for this configuration at which the unstable growth is observed.

The Fig. B.21(a) comprises two y-axes. In the primary y-axis, axial force, *F* has been plotted against the axial displacement, *u* while in secondary y-axis, tangent stiffness Δ*F*/*d*Δ*u* has been depicted against the *u*. In this study, the tangent stiffness is utilised as it represents the stiffness of the laminate at that given instant under the applied compressive load. In Fig. B.21(b) out-of-plane displacements, *w* in the two sublaminates and damage, *d* in cohesive elements have been shown at different axial force levels. In the analysis, the force is ramped up to a reference level of 8 kN. Initially, up to an axial displacement of 0.15 mm, corresponding to an axial force of 5.8 kN, no stiffness reduction is observed. Beyond this point, the stiffness begins to decrease, which can be attributed to geometric nonlinearity caused by the increasing out-of-plane displacement of the thin sublaminate. As shown in Fig. B.21(b), at a force level of 6.720 kN, the thin sublaminate exhibits a displacement of 0.1904 mm, along with some damage at the delamination front. At a force level of 8 kN, a sudden increase in axial displacement occurs due to unstable delamination growth, as shown in Fig. B.21(b). This leads to a significant reduction in stiffness and thus failure of the laminate. Similarly, the axial force and displacement values for each other cases are utilised to calculate the tangent stiffness.
<figure><img src="image_25.png" /><figcaption>Fig. A.20. Schematic representation of the implementation of static CZM (Turon et al. [33]).</figcaption></figure>
<figure><img src="image_26.png" /></figure>

<figure><img src="image_27.png" /><figcaption>Fig. B.21. For near surface delamination configuration with delamination size of 16 mm and imperfection of 0.1%: (a) Primary y-axis shows Axial force versus axial displacement and secondary y-axis shows tangent stiffness drop with respect to axial displacement. (b) The figures highlight out-of-plane displacement, w (left) and damage, d (right) at different force levels. For enhanced visualisation the magnification scale of out-of-plane displacement plots have been increased to 3.0.</figcaption></figure>

## Data availability

The necessary data required to replicate the findings are available in the manuscript.

## References

[1] Bak BLV, Sarrado C, Turon A, Costa J. Delamination under fatigue loads in composite laminates: a review on the observed phenomenology and computational methods. Appl Mech Rev 2014;66(6):060803. http://dx.doi.org/10.1115/1.4027647.

[2] Mekonnen AA, Woo K, Kang M, Kim I-G. Effects of size and location of initial delamination on post-buckling and delamination propagation behavior of laminated composites. Int J Aeronaut Space Sci 2020;21:80-94. http://dx.doi.org/10.1007/s42405-019-00195-0.

[3] Brunner AJ. Experimental study of delamination in cross-ply laminates. In: Delamination behaviour of composites. Woodhead publishing series in composites science and engineering, Elsevier; 2008, p. 281-309. http://dx.doi.org/10.1533/9781845694821.3.281.

[4] Sridharan S. Delamination behaviour of composites. Woodhead Publishing; 2008, p. 1-762.

[5] Mohammadi B, Shahabi F. On computational modeling of postbuckling behavior of composite laminates containing single and multiple through-the-width delaminations using interface elements with cohesive law. Eng Fract Mech 2016;152:88-104. http://dx.doi.org/10.1016/j.engfracmech.2015.04.005.

[6] Pascoe J-A. Slow-growth damage tolerance for fatigue after impact in FRP composites: Why current research won't get us there. Theor Appl Fract Mech 2021;116:103127. http://dx.doi.org/10.1016/j.tafmec.2021.103127.

[7] Whitcomb JD, Shivakumar KN. Strain-energy release rate analysis of plates with postbuckled delaminations. J Compos Mater 1989;23(7):714-34. http://dx.doi.org/10.1177/002199838902300705.

[8] Kyoung W-M, Kim C-G, Hong C-S. Buckling and postbuckling behavior of composite cross-ply laminates with multiple delaminations. Compos Struct 1998;43(4):257-74. http://dx.doi.org/10.1016/S0263-8223(98)00101-9.

[9] Gu H, Chattopadhyay A. An experimental investigation of delamination buckling and postbuckling of composite laminates. Compos Sci Technol 1999;59(6):903-10. http://dx.doi.org/10.1016/S0266-3538(98)00130-4.

[10] Hwang S-F, Liu G-H. Experimental study for buckling and postbuckling behaviors of composite laminates with multiple delaminations. J Reinf Plast Compos 2002;21(4):333-49. http://dx.doi.org/10.1177/0731684402021004253.

[11] Hwang S-F, Huang S-M. Postbuckling behavior of composite laminates with two delaminations under uniaxial compression. Compos Struct 2005;68(2):157-65. http://dx.doi.org/10.1016/j.compstruct.2004.03.010.

[12] Kharazi M, Ovesy HR. Postbuckling behavior of composite plates with through-the-width delaminations. Thin-Walled Struct 2008;46(7):939-46. http://dx.doi.org/10.1016/j.tws.2008.01.005.

[13] Aslan Z, Şahin M. Buckling behavior and compressive failure of composite laminates containing multiple large delaminations. Compos Struct 2009;89(3):382-90. http://dx.doi.org/10.1016/j.compstruct.2008.08.011.

[14] Kharazi M, Ovesy HR, Taghizadeh M. Buckling of the composite laminates containing through-the-width delaminations using different plate theories. Compos Struct 2010;92(5):1176-83. http://dx.doi.org/10.1016/j.compstruct.2009.10.019.
[15] Ovesy HR, Kharazi M. Compressional stability behavior of composite plates with through-the-width and embedded delaminations by using first order shear deformation theory. Comput Struct 2011;89(19):1829-39. http://dx.doi.org/10.1016/j.compstruc.2010.10.016, Civil-Comp.

[16] Park D-Y, Lee S-Y. Parametric effects on embedded delamination buckling in composite structures using the EAS three-dimensional element. Compos Part B: Eng 2012;43(4):1800-12. http://dx.doi.org/10.1016/j.compositesb.2012.02.033.

[17] Gaudenzi P, Perugini P, Riccio A. Post-buckling behavior of composite panels in the presence of unstable delaminations. Compos Struct 2001;51(3):301-9. http://dx.doi.org/10.1016/S0263-8223(00)00146-X.

[18] Wang S, Zhang Y. Buckling, post-buckling and delamination propagation in debonded composite laminates part 2: Numerical applications. Compos Struct 2009;88(1):131-46. http://dx.doi.org/10.1016/j.compstruct.2008.02.012.

[19] Liu PF, Hou SJ, Chu JK, Hu XY, Zhou CL, Liu YL, Zheng JY, Zhao A, Yan L. Finite element analysis of postbuckling and delamination of composite laminates using virtual crack closure technique. Compos Struct 2011;93(6):1549-60. http://dx.doi.org/10.1016/j.compstruct.2010.12.006.

[20] Nikrad SF, Keypoursangsari S, Asadi H, Akbarzadeh AH, Chen ZT. Computational study on compressive instability of composite plates with off-center delaminations. Comput Methods Appl Mech Engrg 2016;310:429-59. http://dx.doi.org/10.1016/j.cma.2016.07.021.

[21] Mekonnen AA, Woo K, Kang M, Kim I-G. Post-buckling and delamination propagation behavior of composite laminates with embedded delamination. J Mech Sci Technol 2020;34(3):1099-110. http://dx.doi.org/10.1007/s12206-020-0211-7.

[22] Andraju LB, Ramji M, Raju G. Snap-buckling and failure studies on CFRP laminate with an embedded circular delamination under flexural loading. Compos Part B: Eng 2021;214:108739. http://dx.doi.org/10.1016/j.compositesb.2021.108739.

[23] Raimondo A, Urcelay Oca I, Bisagni C. Influence of interface ply orientation on delamination growth in composite laminates. J Compos Mater 2021;55(27):3955-72. http://dx.doi.org/10.1177/00219983211031636.

[24] Li Y, Wang B, Zhou L. Study on the effect of delamination defects on the mechanical properties of CFRP composites. Eng Fail Anal 2023;153:107576. http://dx.doi.org/10.1016/j.engfailanal.2023.107576.

[25] Kumar AP, Méndez JP, Burela RG, Bisagni C, Harursampath D, Ponnusami SA. Non-linear buckling analysis of delaminated hat-stringer panels using variational asymptotic method. Compos Struct 2024;345:118276. http://dx.doi.org/10.1016/j.compstruct.2024.118276.

[26] Kutlu Z, Chang F-K. Composite panels containing multiple through-the-width delaminations and subjected to compression. Part II: experiments & verification. Compos Struct 1995;31(4):297-314. http://dx.doi.org/10.1016/0263-8223(95)00093-3.

[27] De Carvalho NV, Krueger R, Mabson GE, Deobald LR. Combining progressive nodal release with the virtual crack closure technique to model fatigue delamination growth without re-meshing. In: 2018 AIAA/aSCE/AHS/aSC structures, structural dynamics, and materials conference. Virginia; 2018.

[28] Dávila CG. From S-N to the Paris law with a new mixed-mode cohesive fatigue model for delamination in composites. Theor Appl Fract Mech 2020;106:102499. http://dx.doi.org/10.1016/j.tafmec.2020.102499.

[29] Akterskaia M, Jansen E, Hallett SR, Weaver P, Rolfes R. Analysis of skin-stringer debonding in composite panels through a two-way global-local method. Compos Struct 2018;202:1280-94. http://dx.doi.org/10.1016/j.compstruct.2018.06.064.

[30] Harper PW, Hallett SR. Cohesive zone length in numerical simulations of composite delamination. Eng Fract Mech 2008;75(16):4774-92. http://dx.doi.org/10.1016/j.engfracmech.2008.06.004.

[31] Donadon MV, de Almeida SFM. Damage modeling in composite structures. In: Comprehensive materials processing. Elsevier; 2014, p. 111-47.

[32] Baluch AH, Falcó O, Jiménez JL, Tijs BHAH, Lopes CS. An efficient numerical approach to the prediction of laminate tolerance to barely visible impact damage. Compos Struct 2019;225:111017. http://dx.doi.org/10.1016/j.compstruct.2019.111017.

[33] Turon A, González EV, Sarrado C, Guillamet G, Maimí P. Accurate simulation of delamination under mixed-mode loading using a cohesive model with a mode-dependent penalty stiffness. Compos Struct 2018;184:506-11. http://dx.doi.org/10.1016/j.compstruct.2017.10.017.

[34] Turon A, Camanho PP, Costa J, Renart J. Accurate simulation of delamination growth under mixed-mode loading using cohesive elements: Definition of interlaminar strengths and elastic stiffness. Compos Struct 2010;92(8):1857-64. http://dx.doi.org/10.1016/j.compstruct.2010.01.012.

[35] SIMULIA User Assistance 2023. ABAQUS/Standard user's manual, version 2023. United States: Dassault Systèmes Simulia Corp; 2023.

[36] Duvant G, John CW, Lions JL. Inequalities in mechanics and physics. Grundlehren der mathematischen Wissenschaften, Springer; 1976.

[37] Camanho PP, Dávila CG, de Moura MF. Numerical simulation of mixed-mode progressive delamination in composite materials. J Compos Mater 2003;37(16):1415-38. http://dx.doi.org/10.1177/0021998303034505.

[38] Martinez X, Oller S, Rastellini F, Barbat AH. A numerical procedure simulating RC structures reinforced with FRP using the serial/parallel mixing theory. Compos Struct 2008;86(15):1604-18. http://dx.doi.org/10.1016/j.compstruc.2008.01.007.

[39] Turon A, Dávila CG, Camanho PP, Costa J. An engineering solution for mesh size effects in the simulation of delamination using cohesive zone models. Eng Fract Mech 2007;74(10):1665-82. http://dx.doi.org/10.1016/j.engfracmech.2006.08.025.

[40] Reeder JR, Demarco K, Whitley KS. The use of doubler reinforcement in delamination toughness testing. Compos Part A: Appl Sci Manuf 2004;35(11):1337-44. http://dx.doi.org/10.1016/j.compositesa.2004.02.021.

[41] Joosten MW, Yang QD, Blacklock M, Wang CH. A cohesive network approach for modelling fibre and matrix damage in composite laminates. Compos Struct 2018;206:658-67. http://dx.doi.org/10.1016/j.compstruct.2018.08.087.

[42] Raimondo A, Doesburg SA, Bisagni C. Numerical study of quasi-static and fatigue delamination growth in a post-buckled composite stiffened panel. Compos Part B: Eng 2020;182:107589. http://dx.doi.org/10.1016/j.compositesb.2019.107589.

[43] Chai H, Babcock CD, Knauss WG. One dimensional modelling of failure in laminated plates by delamination buckling. Int J Solids Struct 1981;17(11):1069-83. http://dx.doi.org/10.1016/0020-7683(81)90014-7.

[44] Köllner A, Kashtalyan M, Guz I, Völlmecke C. On the interaction of delamination buckling and damage growth in cross-ply laminates. Int J Solids Struct 2020;202:912-28. http://dx.doi.org/10.1016/j.ijsolstr.2020.05.035.