# An efficient preselection approach to determine critical parameters for slow delamination growth in postbuckled composites

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

► To cite this version:

S.S. Venkat, P.M. Anilkumar, S. Scheffler, E. Baranger, R. Rolfes. An efficient preselection approach to determine critical parameters for slow delamination growth in postbuckled composites. International Journal of Solids and Structures, 2025, 321, pp.113537. 10.1016/j.ijsolstr.2025.113537. hal-05360321

HAL Id: hal-05360321

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

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
<figure><img src="image_2.png" /></figure>

# An efficient preselection approach to determine critical parameters for slow delamination growth in postbuckled composites

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

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

## ARTICLE INFO

**Keywords:**

- Buckling
- Delamination
- Fatigue
- Imperfection
- Laminates

## ABSTRACT

The presence of delaminations can severely undermine the structural integrity of composite structures. These delaminated sites within a structure are vulnerable to sublimate buckling events due to their reduced bending stiffness. This can in turn progress the delamination growth resulting in loss of load carrying capacity of the structure. The current study presents insights into expected growth behaviour of delaminations in composite laminates under service fatigue loads. Delamination growth can be classified as unstable when the growth rate increases with delamination size or when an abrupt, instantaneous propagation occurs. Conversely, growth is considered stable when the growth rate decreases as the delamination size increases. To gain deeper insights into these growth characteristics, a comprehensive study on composite laminates with through-width delamination has been carried out using an established analytical modelling approach from literature, which is based on the variational principle of total potential energy. For the analysis of delamination growth, the Griffith criterion is used, which postulates that energy release rate (G) governs the process. It has been shown that G begins to increase following thin sublimate buckling, with the rate of increase with respect to applied strain being higher for a laminate with short delamination than for one with long delamination. The expected delamination growth characteristics under service fatigue loading has been analysed, and strain allowables corresponding to the fatigue limits have been determined. Similarly for static loading scenarios, conditions for stable and unstable growth have been presented and linked to fatigue growth characteristics. A critical discussion is subsequently presented on the effects of imperfection magnitude and delamination depth on the growth behaviour. By presenting these insights into growth behaviour, this study proposes an efficient approach for identifying key parameters that would ensure slow growth of delaminations in postbuckled composite laminates.

## 1. Introduction

Owing to their exceptional specific stiffness and strength, composite materials have become a preferred choice for a wide range of structural applications, spanning the aerospace and wind energy industries. However, the operational reliability of composite structures can be significantly compromised by various damage mechanisms, with delaminations being one of the major concern (Sridharan, 2008). Delaminations can be caused by impact events or be triggered by static or fatigue loads, often due to matrix cracking, stress concentrations at material or geometric discontinuities, or interfacial defects introduced during manufacturing (Sridharan, 2008; Brunner, 2008). Composite structures with these delaminations are highly prone to premature failure, especially under compressive loads (Mohammadi and Shahabi, 2016; Mekonnen et al., 2020a). This occurs because delaminations locally reduce the bending stiffness of the structures, making them susceptible to buckling and subsequent growth. While the buckling event

itself can lead to a loss of structural stiffness, the subsequent growth of delaminations can drastically degrade the stiffness and consequently, the load-carrying capacity of the structure. As a result, this problem has attracted considerable attention over the years, with several numerical parametric studies conducted to examine the influence of factors such as delamination size, location and shape (Nikrad et al., 2016; Mekonnen et al., 2020b). One of the key findings of these parametric studies is that composite laminates with a single through-width delamination can experience three distinct buckling modes as shown in Fig. 1. Here, a through-width delamination refers to a delamination in a composite column or strut, where the delamination extends across the entire width of the laminate.

When delamination is located near the laminate surface, known as a near-surface configuration, it locally divides the laminate into thin and thick sublaminates. For this case the laminate can undergo either local

<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.
E-mail address: s.venkat@isd.uni-hannover.de (S.S. Venkat).</footnote>
<figure><img src="image_4.png" /><figcaption>**Fig. 1.** Buckling modes in composite laminates with a through-width delamination (Mekonnen et al., 2020a; Venkat et al., 2025).</figcaption></figure>

or a mixed local/global buckling mode (Mekonnen et al., 2020a). The local buckling mode occurs when the stiffness of thick sublamine is significantly greater than the stiffness of the thin sublamine, causing the thin sublamine to undergo buckling while the thick sublamine undergoes negligible bending in opposite direction. The second buckling mode is mixed local/global which occurs when the ratio of stiffness of thick sublamine to thin sublamine is considerably lower than in local buckling case. In this mode, in addition to thin sublamine buckling, the thick sublamine also exhibits bending in the opposite direction. In case of mid-plane delamination where both the sublaminates generally have similar stiffness, the laminate tends to undergo global buckling, similar to what is expected in an intact laminate (Venkat et al., 2025; Mekonnen et al., 2020a). In a given laminate, local and mixed local/global modes are more critical, as they can occur at much lower applied strains than the global mode. This is due to the significantly lower bending stiffness of the thin sublamine. In another study, Köllner et al. (2020) have presented a detailed analysis of the dependence of buckling responses on the length and depth of through-width delaminations in thin cross-ply laminates, using a Rayleigh-Ritz-based analytical modelling approach. It has been reported that, for very short delaminations, the laminate exhibits a global closed buckling mode after reaching a limit point. With a slight increase in delamination length, the laminate initially displays a global closed-mode response, which transitions to an opening mode after reaching a limit point. However, for longer delaminations, the laminates exhibit opening mode with subsequent global response or only opening mode throughout the deformation history. It has also been shown that the delamination lengths corresponding to these responses vary with the depth of the delamination.

However, as mentioned before, it is crucial to not only understand the buckling modes of composite laminates with delaminations, but also to analyse their subsequent growth, especially the growth characteristics. In an early study, Chai et al. (1981) examined the quasi-static delamination growth characteristics of composite laminates with predefined delamination. For quasi-static loading the growth is stable when the applied strains have to be increased for the growth to progress and it is unstable when the delamination grows abruptly to a larger size at the same strain level. However, during their operational lifetime, composite structures are subjected to service loads, such as fatigue, which can cause the growth of delaminations (Bisagni et al., 2013; Bak et al., 2014) and therefore, it is also essential for designers to understand the growth characteristics under service fatigue loads. For the certification requirements, there are two well known damage tolerance based design concepts for the composite structures (Pascoe, 2021). In the more prevalent ‘no-growth’ design concept the strain allowables are kept below a low value known as fatigue threshold,

ensuring that existing damages, such as delaminations, does not propagate (Pascoe, 2020). This approach is however highly conservative and leads to a substantial weight penalty. The ‘slow-growth’ design concept on the other hand allows the existing damages originating from impact or manufacturing imperfections to grow, as long as the growth remains slow, stable and predictable. This approach is less conservative compared to ‘no-growth’ approach as the strain allowable can be increased beyond the fatigue threshold as long the growth remains stable and slow. Stable growth requires that, under a constant fatigue load, the growth rate does not increase as the delamination size increases (Chai et al., 1981; Köllner et al., 2020). However, at the same time, it must also be ensured that the residual stiffness of the structure is not compromised with the damage growth, hindering its effective functioning. The requirement for damage growth to be slow and stable is to effectively define the inspection window. It must also be predictable, allowing designers to setup an effective inspection schedule with adequate checks in that window. In order to design laminates for delamination buckling using a less conservative ‘slow-growth’ approach, initially a comprehensive parametric study of the growth characteristics is a requirement. These growth characteristics under quasi-static and fatigue loading share certain similarities, particularly in the one of the definition of unstable growth as the abrupt extension of delamination under applied constant strain. Therefore, an initial comprehensive parametric study under quasi-static loading can also offer valuable insights which are applicable to fatigue loading. Such a quasi-static parametric study would involve various factors like delamination size, sublamine stiffness ratios, applied strain levels and the geometric imperfection magnitudes.

One possibility of performing this would be using numerical methods from literature. The primary approach in literature for predicting buckling point and postbuckling deformation has been the geometrically nonlinear finite element (FE) method (Liu et al., 2011; Mohammadi and Shahabi, 2016; Nikrad et al., 2016; Mekonnen et al., 2020a,b; Venkat et al., 2025). It has often been combined with methods for delamination propagation such as fracture mechanics based virtual crack closure technique (VCCT) (Liu et al., 2011) or damage mechanics-based cohesive zone model (CZM) (Mohammadi and Shahabi, 2016; Nikrad et al., 2016; Mekonnen et al., 2020a,b; Venkat et al., 2025). Recently, Venkat et al. (2025) conducted a parametric study employing CZM within a geometrically nonlinear FE framework. A quasi-static ‘Effect of Defect’ study was performed on composite laminates to evaluate whether a pre-existing through-width delamination defect could exhibit unstable abrupt growth under a constant load. The influence of through-thickness location was investigated by maintaining a constant load while systematically increasing the length of delamination in the FE models. Other studies in literature have mainly focused on determining the critical failure load in composite laminates associated with various delamination related factors. These include through-thickness (Mohammadi and Shahabi, 2016; Mekonnen et al., 2020a) and axial locations (Nikrad et al., 2016), delamination length (Mekonnen et al., 2020a; Mohammadi and Shahabi, 2016) and different delamination types (Wang and Zhang, 2009; Liu et al., 2011; Mohammadi and Shahabi, 2016; Mekonnen et al., 2020a,b). The delamination types include through-width delamination (Wang and Zhang, 2009; Liu et al., 2011; Mohammadi and Shahabi, 2016; Mekonnen et al., 2020a), which has been previously defined, and another type called embedded delamination (Mekonnen et al., 2020b) in composite plates. In embedded type, delaminations of various shapes, such as rectangle or circle, are contained within the laminate and do not span its full width. Studies have also been performed to analyse the buckling and postbuckling delamination growth in co-cured composite stiffened panels (Zhang et al., 2024; Raimondo et al., 2020) and unstiffened curved composite panels (Ameri et al., 2020). Although initial studies can be performed using numerical methods, their advanced model setup effort and computational costs, make them inefficient for a comprehensive study, such as the one planned in this work. The
computational costs are mainly due to the need for a fine mesh for delamination propagation methods, along with convergence difficulties due to stiffness degradation in implicit analysis, that reduce the load increment size or increase the number of Newton-Raphson iterations required for convergence at each increment. As a solution to this challenge, efficient analytical modelling approaches discussed hereafter can provide valuable insights within a short time, thereby significantly reducing the effort compared to numerical analyses.

Analytical methods to gain insights into quasi-static growth behaviour following sublamine buckling in composites have been explored since the early 1980s. This research began by the seminal work of Chai et al. (1981), which focused on through-width near-surface delamination configurations. In this work initially a so-called 'Thin Film' model was developed to address a theoretical laminate in which the thick sublamine does not undergo bending. This was done by considering the thick sublamine of infinite thickness which greatly simplified the formulation. A more cumbersome model based on beam theory was also presented for laminates of finite thickness, which incorporated the effects of bending in the thick sublamine. In this work the laminate was considered to be linear elastic and isotropic in both the models. The delamination growth criterion was governed by fracture mechanics-based energy release rate G. In a later work by Chai and Babcock (1985) an analytical method was developed to study growth criterion of embedded elliptical near surface delamination in a laminate. This was done similarly by breaking the overall problem into an elastic stability problem and a fracture problem based on G. Later on studies were performed by Sheinman and Soffer (1991), and Sheinman et al. (1998) on delamination buckling. In their first study (Sheinman and Soffer, 1991) postbuckling response of a through-width delaminated beam was studied without considering a growth criterion. Nonlinear analysis based on Von Karman kinematic approach was considered to study the influence of different boundary conditions (simply supported and clamped), stacking sequences and imperfections. Similarly in the later study by Sheinman et al. (1998), parametric studies were performed by using variational principle to derive the equilibrium equations for the problem and by including a growth criterion based on G. Other early studies have focused exclusively on determining buckling loads in delaminated composite laminates, without addressing the postbuckling deformation response or subsequent delamination growth criteria (Simitses et al., 1985; Huang and Kardomateas, 1998).

Furthermore in several studies (Kharazi and Ovesy, 2008; Hunt et al., 2004), researchers developed various analytical methods for predicting the buckling load and then postbuckling deformation in composite laminates with delaminations. In an initial study by Kharazi and Ovesy (2008), they applied classical laminated plate theory (CLPT) and utilised a Rayleigh-Ritz based formulation to investigate single near-surface delamination. Building on this, their subsequent work compared CLPT with first-order and higher-order shear deformation theories (FSDT and HSDT) (Kharazi et al., 2010). In a specific study by Ovesy and Kharazi (2011), they explored the effects of two through-width delaminations and a single circular embedded delamination using FSDT. However, their studies were limited to deformation behaviour and a growth criterion was not defined. In recent years Köllner et al. proposed advances in analytical modelling approaches in laminates with a single through-width delamination in a composite strut (Köllner et al., 2016, 2020) and a single elliptical delamination in a composite plate (Köllner and Völlmecke, 2018). In their studies, the pre- and postbuckling deformation paths were determined based on the well known variational principle of total potential energy (TPE) and then the fracture criterion based on G was used. For through-width delamination, the proposed formulation extended beyond traditional fracture criterion by also accounting for deformation behaviour after the delamination growth was triggered under quasi-static loading conditions (Köllner et al., 2016). This was achieved through the use of an extended total potential energy (ETPE) formulation. In a later study by Köllner et al. (2020), matrix crack density was accounted for in

a cross-ply laminate by considering an equivalent constrained model (ECM) within the delamination buckling problem. The ECM was used to determine reduced stiffness properties based on the matrix crack densities. In summary, efficient analytical methods have been extensively proposed in the literature to address various aspects of the delamination buckling problem, including buckling load, postbuckling response, quasi-static growth criteria, and deformation behaviour after delamination growth under quasi-static loading conditions. These methods have also been applied to different types of delaminations, such as through-width and embedded delaminations. Moreover, delamination growth characteristics have also been examined in earlier studies, such as the work by Chai et al. (1981), in which investigations were done for local and mixed local/global buckling modes in isotropic delaminated beams. More recently, Köllner et al. (2020) examined cross-ply laminates, linking delamination growth characteristics to the previously discussed distinct buckling responses identified in their study. However, to gain a more generalised understanding of delamination growth characteristics in postbuckled laminates, it is also essential to investigate other laminate types, such as quasi-isotropic configurations, which are the focus of this study. Furthermore, as this work is primarily concerned with fatigue loading, the main objective is to identify the key parameters that govern slow delamination growth behaviour under such loading conditions. Therefore, a comprehensive analysis is needed to investigate how various parameters such as loading, geometric imperfections, sublamine stiffness ratios and delamination depth affect fatigue growth characteristics.

In the present study, the comprehensive investigations have been carried out in laminates with a through-width delamination. The analysis considers various above mentioned factors, including applied strain levels, imperfection magnitude and sublamine stiffness ratios incorporating also the influence of delamination depth. The novelty of this study lies in its focus on developing a comprehensive understanding of the expected delamination growth behaviour in a quasi-isotropic laminate subjected to fatigue loading conditions. These insights will serve as a crucial step towards defining less conservative fatigue strain allowables for the problem of delamination buckling and identifying the conditions that enable slow and stable delamination growth, thereby offering an alternative to the traditional 'no-growth' damage tolerance approach. To conduct this parametric study, the analytical modelling approach developed by Köllner et al. (2020) has been utilised. This modelling approach incorporates the ability to account for arbitrary stacking sequences of composite laminates, define geometric imperfections, and apply the Griffith failure criterion based on G.

The paper is organised into four sections. Section 2 outlines the methodology, including a description of the analytical modelling framework, and an example comparing the analytical results with FE to verify the implementation. Section 3 starts with a general aspects of the conducted study followed by discussion of the results. The paper ends with Section 4, which offers the conclusions and an outlook to future work.

## 2. Methodology

### 2.1. Modelling approach

This subsection introduces an analytical model used to investigate the effects of delamination defects on postbuckling response and the associated growth characteristics. The model, based on the work of Köllner et al. (2020), employs the variational principle of the total potential energy to evaluate the postbuckling response. The out-of-plane displacements in the model are approximated using Rayleigh-Ritz method and it additionally incorporates a fracture criterion. The flowchart below (Fig. 2(a)) outlines the steps in implementation of the analytical modelling framework. On the right side of the flowchart, corresponding illustrations (Fig. 2(b)–(e)) are provided for each step. To analyse the buckling and postbuckling response of the delaminated
composite laminate, the laminate is subdivided into three parts as shown. Parts (1) and (2) are the thin and thick sublaminates in the delaminated region, while the other part is identical intact region of the laminate located axially on either side of the delamination, which is denoted as part (3). Hereafter, the superscript (k) in the variables represents each of the three parts of the laminate (k = 1 to 3). Initially all inputs to the analytical model (material and geometric parameters, delamination length (l) and end shortening displacement (Δ)) are to be defined, including specifying the imperfection magnitude (q₀) for the thin sublaminate in the delaminated region. In this formulation, the imperfection is represented as a sine function, ensuring that at the two ends corresponding to the delamination front, both out-of-plane displacements and slopes are zero. For the sake of brevity, only the key aspects of the formulation are highlighted here, while further important formulations are comprehended in Appendix. It is recommended to refer to Köllner et al. (2020) for in-depth details on the adopted model.

The postbuckled deformation of the laminate is described by the following four unknown generalised coordinates in this formulation:

* q₁: The out-of-plane deformation amplitude of part (1)

* q₂: The out-of-plane deformation amplitude of part (2)

* q₃: Slope or rotation at the interface between intact and delaminated regions

* q₄: The end-shortening of the delaminated region accounts for contributions from both axial shortening and out-of-plane deformation.

The out-of-plane buckling displacements ($w^{(k)}$) in each laminate part is a function of the above mentioned generalised coordinates. For example it is evident that the buckling displacement of part (2) is influenced not only by $q_2$ but also by $q_3$. In this model the buckling displacements have been approximated using a Rayleigh-Ritz formulation considering the geometric boundary and continuity conditions. These conditions along with the displacement equations have also been presented in Appendix.

The CLPT has been applied to account for arbitrary composite layups. To account for the various coupling terms of the ABD matrix derived from CLPT, effective stiffness parameters $A_{\text{eff}}^{(k)}$ for in-plane, $B_{\text{eff}}^{(k)}$ for coupling, and $D_{\text{eff}}^{(k)}$ for bending are defined for different parts of the laminate, as described in Köllner et al. (2016). To determine effective stiffnesses, initially the forces $n_{yy}$ and $n_{xy}$ and moments $m_{yy}$ and $m_{xy}$ in the CLPT formulation are assumed to be zero. For this problem of a laminate strut with through-width delamination, the assumption is reasonable as mainly the force $n_{xx}$ and moment $m_{xx}$ acts on various parts of the laminate. By taking this assumption, the resulting strains ($\epsilon_{yy}, \epsilon_{xy}$) and curvatures ($\kappa_{yy}, \kappa_{xy}$) are determined in terms of elements of ABD matrix and $\epsilon_{xx}$ and $\kappa_{xx}$. These are then substituted in the expressions for $n_{xx}$ and $m_{xx}$ which results in effective stiffness terms. These effective stiffnesses thereby include the effects of various coupling terms: stretching-shearing ($A_{16}, A_{26}$), stretching-stretching ($A_{12}$), stretching-twisting ($B_{16}, B_{26}$), stretching-bending ($B_{12}$), bending-twisting ($D_{16}, D_{26}$) and bending-bending ($D_{12}$) interactions. The in-plane strains in all the parts of laminate ($\epsilon_0^{(k)}$) are defined using a non-linear kinematic approach based on Hunt et al. (2004). These in-plane strains are determined by dividing the corresponding axial shortening terms ($u^{(k)}$) by the respective length of the parts. The expressions for axial shortening terms have also been summarised in Appendix.

The curvatures for each part is given by $\kappa_0^{(k)} = \frac{\partial^2 w^{(k)}}{\partial x_k^2}$. The model employs variational principle of TPE to determine the out-of-plane displacements with respect to applied end-shortening displacement (Δ). For this, the total strain energy (U) is computed as sum of contributions from the membrane ($U_{\text{Mem}}^{(k)}$), bending ($U_{\text{Bend}}^{(k)}$) and coupled membrane-bending ($U_{\text{Coupled}}^{(k)}$) terms for each of the three parts of the laminate, as shown below,

$$U = \sum_{k=1}^{3} \left( U_{\text{Mem}}^{(k)} + U_{\text{Bend}}^{(k)} + U_{\text{Coupled}}^{(k)} \right) \quad (1)$$

where the individual energy contributions for the first two parts of the laminate (k = 1 to 2) is given by,

$$ \begin{aligned} U_{\text{Mem}}^{(k)} &= \frac{1}{2} b \int_{0}^{l} A_{\text{eff}}^{(k)} (\epsilon_{0}^{(k)})^{2} dx_{k}, \\ U_{\text{Bend}}^{(k)} &= \frac{1}{2} b \int_{0}^{l} D_{\text{eff}}^{(k)} (\kappa_{0}^{(k)})^{2} dx_{k}, \\ U_{\text{Coupled}}^{(k)} &= \frac{1}{2} b \int_{0}^{l} 2 B_{\text{eff}}^{(k)} \kappa_{0}^{(k)} \epsilon_{0}^{(k)} dx_{k} \end{aligned} \quad (2) $$

For the part (3), the energies for only one intact region is calculated by applying integration limits from 0 to $l_1 = sl$ and then doubled (for left and right part) as shown below:

$$ \begin{aligned} U_{\text{Mem}}^{(3)} &= b \int_{0}^{sl} A_{\text{eff}}^{(3)} (\epsilon_{0}^{(3)})^{2} dx_{3}, \\ U_{\text{Bend}}^{(3)} &= b \int_{0}^{sl} D_{\text{eff}}^{(3)} (\kappa_{0}^{(3)})^{2} dx_{3}, \\ U_{\text{Coupled}}^{(3)} &= b \int_{0}^{sl} 2 B_{\text{eff}}^{(3)} \kappa_{0}^{(3)} \epsilon_{0}^{(3)} dx_{3} \end{aligned} \quad (3) $$

Under end-shortening displacement loading, the TPE is equal to the total strain energy, $\Pi = U$ (Kharazi and Ovesy, 2008). At equilibrium, the first variation of the TPE equals zero ($\delta\Pi = 0$), leading to the condition $\frac{\partial H}{\partial q_i} = 0.0$ for $i = 1$ to 4. This results in a system of nonlinear equations, which is solved numerically using the Newton-Raphson algorithm to determine the unknowns, $q_i$. The calculation has been carried out using the software, MATLAB. Once these unknowns are obtained, the buckling displacement $w^{(k)}(q_i)$ for each part can be computed.

For the analysis of delamination growth, a Griffith-type quasi-brittle fracture process is assumed, which is based on G. G is determined using the equation shown in the flowchart (Fig. 2). Under quasi-static loading, delamination growth occurs when $G \ge G_c$. Mode I $G_c$ is adopted in this study as a conservative assumption, given that for graphite/epoxy composite materials, Mode I $G_c$ is typically lower than Mode II $G_c$ (Harper and Hallett, 2008). Under fatigue loading, the growth takes place when, $G_{\text{th}} < G < G_c$, where $G_{\text{th}}$ represents the lower threshold, below which no growth occurs (Bak et al., 2014). The dependence of crack growth rate on the G is governed by the general form of Paris law as shown in Fig. 2 (Bak et al., 2014). In this study $f(G) = G_{\text{max}}$ is used where $G_{\text{max}}$ represents the maximum energy release rate. The influence of the minimum energy release rate, $G_{\text{min}}$ (or load ratio) is implicitly incorporated through the parameters C and m, which also account for other factors such as the fracture mode mixity.

## 2.2. Verification example

The analytical model adopted from Köllner et al. (2020) (see subsection 2.1) has been implemented, and to confirm the accuracy of the implementation, the results are first verified against the FE results reported in the same reference. In this example a slender composite laminate column with thickness of 2.49 mm and a layup of [(0/90)$_{s}$]$_{l7}$ is chosen. The postbuckling deformation response and the critical force associated with quasi-static delamination growth in the selected laminate with a single through-width delamination have been analysed. The laminate has a total of 28 plies with a delamination placed in between plies 23 and 24. The results shown in Fig. 3 correspond to three different initial lengths of delamination. Details of geometric and material parameters of the laminates are taken from Köllner et al. (2020). The FE predictions used for this verification have been conducted in Köllner et al. (2020) using ABAQUS software. The plies have been modelled using layered shell elements and the VCCT has been employed to calculate G. Similar to analytical modelling approach, the FE simulations also utilised a mode I crack growth criterion where the growth occurs when $G \ge G_c$. In Fig. 3, the delamination length is normalised with respect to laminate length ($l_{\text{norm}}$), the force is normalised with respect
<figure><img src="image_5.png" /><figcaption>Fig. 2. Schematic representation of analytical framework used in this study.</figcaption></figure>

to the critical buckling force of the intact laminate as ($P_{\text{norm}}$), and the out-of-plane buckling displacements are normalised with respect to the laminate thickness ($w_{\text{norm}}$). Since the analytical model is formulated for displacement based loading, forces are determined by computing the reaction forces in the intact part 3 (Subsection 2.1) of the laminate. The graphs ($P_{\text{norm}}$) versus ($w_{\text{norm}}$) are plotted for both FE and analytical model until failure when $G = G_{\text{IC}}$.

The results corresponding to $l_{\text{norm}} = 0.64$ is explained first. In this configuration, $w_{\text{norm}}$ for the thin sublamine, represented by a dashed blue line, starts to increase when $P_{\text{norm}}$ reaches around 0.1. At this stage, the thin sublamine is considered to have undergone buckling. Similarly, $w_{\text{norm}}$ for the thick sublamine, depicted by a solid blue line, begins to increase significantly once $P_{\text{norm}}$ exceeds 0.2. Due to the slender nature of the laminate, it exhibits a mixed local/global buckling mode. The $w_{\text{norm}}$ continues to increase in thin and thick sublaminates until reaching values of 0.819 and -0.766 respectively, at which point the delamination growth criterion is reached. The FE results, represented by blue markers, show very close overall agreement with the analytical results, including the last marker point which is the failure point. The analytical results for $l_{\text{norm}} = 0.52$ and $l_{\text{norm}} =$

<figure><img src="image_6.png" /><figcaption>Fig. 3. Normalised out-of-plane displacements in thin and thick sublaminates with respect to normalised force. Verification of implementation of the analytical model with the FE results of Köllner et al. (2020).</figcaption></figure>
<table><caption>**Table 1**
Details of laminate dimensions, layup, delamination details adopted in this study.</caption><thead><tr><th rowspan="2">Identification</th><th colspan="2">Details</th></tr><tr><th>Thick laminate</th><th>Thin laminate</th></tr></thead><tbody><tr><td>Laminate layup</td><td>[0<sub>3</sub>/ - 45<sub>3</sub>/45<sub>3</sub>/90<sub>3</sub>]<sub>6s</sub></td><td>[0<sub>3</sub>/ - 45<sub>3</sub>/45<sub>3</sub>/90<sub>3</sub>]<sub>2s</sub></td></tr><tr><td>Laminate thickness (t)</td><td>18.864 mm</td><td>6.288 mm</td></tr><tr><td>Laminate length (L)</td><td colspan="2">100 mm</td></tr><tr><td>Laminate width (b)</td><td colspan="2">10 mm</td></tr><tr><td>Delamination length (l)</td><td colspan="2">Starting with 10 mm</td></tr><tr><td>Delamination location</td><td colspan="2">Between final 0<sub>3</sub> and penultimate -45<sub>3</sub> layers</td></tr><tr><td>End shortening displacement (Δ)</td><td colspan="2">Study on influence</td></tr><tr><td>Imperfection magnitude (q<sub>0</sub>)</td><td colspan="2">0.001, 0.01 and 0.02</td></tr></tbody></table>

<table><caption>**Table 2**
Linear-elastic ply properties and interface fracture properties of IM7/8552 unidirectional graphite/epoxy composite (Joosten et al., 2018).</caption><caption>Ply elastic and interface critical energy release rate properties</caption><thead><tr><th>E<sub>11</sub> [MPa]</th><th>E<sub>22</sub> [MPa]</th><th>ν<sub>12</sub> [-]</th><th>G<sub>12</sub> [MPa]</th><th>G<sub>IC</sub> [N/mm]</th><th>G<sub>ICC</sub> [N/mm]</th></tr></thead><tbody><tr><td>160580</td><td>11 500</td><td>0.31</td><td>5430</td><td>0.21</td><td>0.77</td></tr></tbody></table>

<figure><img src="image_7.png" /><figcaption>**Fig. 4.** Schematic representation of the laminate featuring a near-surface and through-width central delamination.</figcaption></figure>

0.33 similarly demonstrate very good agreement with the FE results. In these cases, the short delamination length delays the buckling of the thin sublaminate and, consequently, also the thick sublaminate buckling. The similar out-of-plane displacements in the thick and thin sublaminates result because of their relatively small thickness ratio of approximately 4.5. In the following sections, a more pronounced contrast will be observed, as few of the evaluated configurations exhibit thickness ratios of 15 or higher. In summary, the results of the analytical model demonstrate the accuracy of its implementation through very good agreement with the FE results, thereby confirming its suitability for the parametric studies planned in this work. While previous studies in literature (Chai et al., 1981; Köllner et al., 2020) have investigated delamination growth behaviour in simplified isotropic beam models or cross-ply laminates, the present study focuses on a quasi-isotropic laminate, with particular emphasis on insights related to fatigue limits.

# 3. Results and discussions

## 3.1. General aspects of the study

Two laminates with different thicknesses, each featuring a near-surface, through-width central delamination, have been studied. Fig. 4 provides a schematic representation of the laminate, with corresponding details such as dimensions, layup configuration, and loading listed in Table 1.

The study focuses on generic examples of quasi-isotropic (QI) laminates: a thick laminate consisting of 48 layers and a thin laminate comprising 16 layers. Each layer of the laminate is made up of three plies (e.g.: $0_3$), and the through-width delamination is located between the penultimate and final layers in each laminate, effectively dividing them into thick and thin sublaminates. As a result, the ratio of thickness of thick-to-thin sublaminates are 47:1 and 15:1 for the thick and thin laminates, respectively. The layups have been chosen for their simplicity, with the sub-layup [0<sub>3</sub>, -45<sub>3</sub>, 45<sub>3</sub>, 90<sub>3</sub>] serving as the basic repeating

unit in the overall laminate structure. In the thin laminate, this unit repeats four times, while in the thick laminate, it repeats twelve times. In both the laminate configurations, the repeating units are arranged symmetrically about the laminate midplane. Three different imperfection magnitudes ($q_0$) have been analysed (Table 1), which represents maximum out-of-plane geometric perturbation at the centre of thin sublaminates. Since strain based loading is considered, the end-shortening strain ($\epsilon_0$) has been defined as the applied displacement ($\Delta$) divided by the initial length of the laminate (L). Various end-shortening strains have been investigated, which are crucial for determining the strain allowable associated with fatigue limits and identifying the conditions that enable stable delamination growth under such loading conditions. Furthermore, in the case of the thin laminate, a critical discussion is presented on how the through-thickness location of the delamination influences the buckling response and growth characteristics. In addition to the first configuration, where the delamination is located between layers 15 and 16 (see Table 2), three additional configurations are examined with the delamination located between layers 14 and 15, 13 and 14, and finally between layers 11 and 12. The study has been performed on a unidirectional graphite/epoxy composite (IM7/8552). For the selected composite, the linear-elastic, transversely isotropic ply properties and the interface critical energy release rates corresponding to fracture modes I and II are given in Table 2.

Various insights into the expected behaviour of the selected composite laminates under compressive fatigue loading are presented in the following subsections. It begins with an evaluation of the sensitivity of delamination growth to loading, followed by an examination of growth characteristics under both fatigue and quasi-static loading, with a link drawn between the two. Subsequently, the influence of imperfection magnitude is examined in the case of thick laminate, while the influence of through thickness location of delamination is analysed in the case of thin laminate. The study concludes with findings from quasi-static FE analysis, providing deeper insights into the advantages of the adopted analytical method and outlining potential future directions for this research.

## 3.2. Discussions on thick laminate

The thick laminate used in this study refers to a laminate with the layup of [0<sub>3</sub>/ - 45<sub>3</sub>/45<sub>3</sub>/90<sub>3</sub>]<sub>6s</sub>. All results discussed below are for this thick laminate with near-surface through-width delamination, located between the final $0_3$ layer and the penultimate -45<sub>3</sub> layer. In the first two subsections, which address the sensitivity of growth to loading and the growth characteristics, an imperfection magnitude ($q_0$) of 0.001 is considered. The final subsection deepens the understanding of growth characteristics by examining the influence of imperfection magnitudes.
### 3.2.1. Insight into sensitivity of delamination growth to load

In this subsection, the insights into sensitivity of delamination growth to fatigue loading have been presented. This has been accomplished by examining the rate of increase of $G$ with applied end-shortening strain ($\epsilon_0$) under quasi-static loading conditions. To study the influence of different delamination lengths, a comparison is made between lengths of 10 mm and 26 mm. The results pertaining to delamination length of 10 mm are presented in Fig. 5. In Fig. 5(a), two y-axes are displayed: the primary y-axis on the left shows the energies ($U$), while the secondary y-axis on the right represents the normalised energy release rate ($G/G_c$) and they are both plotted against $\epsilon_0$ on the x-axis. Energies corresponding to part (1), the thin sublamine in delaminated region have been plotted. They comprise membrane ($U_{Mem}^{(1)}$) and bending ($U_{Bend}^{(1)}$) contributions. Fig. 5(b) shows the graph of out-of-plane displacements at the centres of thin sublamine ($w^{(1)}$) and thick sublamine ($w^{(2)}$), plotted against ($\epsilon_0$).

As shown in Fig. 5(a), when $\epsilon_0$ is initially increased, only $U_{Mem}^{(1)}$ increases, while $U_{Bend}^{(1)}$ remains zero. This indicates, that at this stage, part (1) is undergoing only axial shortening under compression. It is also observed that $G/G_c$ remains zero at this stage. This axial shortening of part (1) continues until $\epsilon_0$ reaches a critical value, $\epsilon_{cr}$ which represents the buckling strain of part (1). Beyond $\epsilon_{cr}$, as expected, $w^{(1)}$ begins to increase, as observed in Fig. 5(b). This increase in $w^{(1)}$ leads to increase in bending contribution $U_{Bend}^{(1)}$ in part (1) while $U_{Mem}^{(1)}$ starts to plateau. It is also observed in Fig. 5(a) that for $\epsilon_0 > \epsilon_{cr}$, $G/G_c$ begins to increase, which is due to the out-of-plane deformation of part (1). Subsequently, the failure strain $\epsilon_f$ under quasi-static loading condition is reached when $G/G_c = 1.0$. For this configuration, $\epsilon_{cr} = 0.0050$ while $\epsilon_f = 0.00543$, which shows that for 8.6% increase in $\epsilon_0$, $G/G_c$ increases from zero to one. At $\epsilon_f$, the out-of-plane deformations $w^{(1)}$ equals $0.123$ mm while $w^{(2)}$ is negligible due to substantial thickness of thick sublamine. This indicates that the configuration exhibits local buckling. To examine how the behaviour varies with delamination length, a long delamination of 26 mm is analysed next as an illustrative example.

For laminate with initial delamination length of 26 mm, the behaviour (Fig. 6) has many similarities to the one with 10 mm (Fig. 5). Initially, due to axial shortening, $U_{Mem}^{(1)}$ increases with $\epsilon_0$ until the $\epsilon_0 = \epsilon_{cr}$. Beyond this point, $U_{Bend}^{(1)}$ begins to increase due to increase in out-of-plane deformation of part (1) ($w^{(1)}$), while $U_{Mem}^{(1)}$ shows minimal increase with $\epsilon_0$. In these aspects, the behaviour is similar to the laminate with a 10 mm delamination length, however, there are also notable differences, which are discussed next. In this case, $\epsilon_{cr}$ is considerably lower, at 0.0007, compared to the short delamination length of 10 mm, where it is 0.0050. The $\epsilon_f$ for this configuration is 0.00219, which is 212.86% higher than its $\epsilon_{cr}$. This represents a significant increase compared to the 10 mm delamination length, where the difference is only 8.6%. Therefore, for long delaminations, the rate of increase of $G$ is considerably more gradual with $\epsilon_0$, as compared to short delamination length of 10 mm. Other major difference is that at $\epsilon_f$, $w^{(1)}$ has a much higher value of 0.650 mm while $w^{(2)}$ is still negligible.

To gain insights into fatigue loading, it is first essential to establish the conditions under which growth occurs during fatigue. The fatigue loading is characterised by a maximum strain (i.e. maximum load) and a load ratio which defines the minimum strain. As noted in subsection 2.1, the form of Paris law adopted in this study is $\frac{dI}{dN} = CG_{max}^m$, where the delamination growth rate is dependent on $G_{max}$ and parameters $C$ and $m$ account for effects such as load ratio and fracture mode mixity. Consequently, this formulation accounts for both the maximum strain through $G_{max}$ and for the load ratio through model parameters. In the context of fatigue related discussions, $\epsilon_0$ in figures above represents maximum end-shortening strain and $G$ represents $G_{max}$. This interpretation also applies to the subsequent sections of this document. It is also mentioned in subsection 2.1, that growth under fatigue occurs for condition $G_{th} < G < G_c$ where $G_{th}$ represents the very low threshold

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

<figure><img src="image_9.png" /><figcaption>Fig. 5. For thick laminate and delamination length of 10 mm: (a) Membrane ($U_{Mem}^{(1)}$) and bending ($U_{Bend}^{(1)}$) energies in the thin sublamine and normalised energy release rate ($G/G_c$) with respect to end shortening strain ($\epsilon_0$). (b) Out-of-plane displacements of thin ($w^{(1)}$) and thick ($w^{(2)}$) sublaminates with respect to $\epsilon_0$.</figcaption></figure>

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

<figure><img src="image_11.png" /><figcaption>Fig. 6. For thick laminate and delamination length of 26 mm: (a) Membrane ($U_{Mem}^{(1)}$) and bending ($U_{Bend}^{(1)}$) energies in the thin sublamine and normalised energy release rate ($G/G_c$) with respect to end shortening strain ($\epsilon_0$). (b) Out-of-plane displacements of thin ($w^{(1)}$) and thick ($w^{(2)}$) sublaminates with respect to $\epsilon_0$.</figcaption></figure>
value of G, below which no growth takes place. This implies that for the growth to occur under fatigue, $\epsilon_0$ must surpass a strain value higher than $\epsilon_{cr}$, ensuring that G is higher than $G_{th}$. However, it must remain below $\epsilon_f$, when $G = G_c$, which corresponds to the quasi-static growth condition. As demonstrated in this subsection, the rate of increase of G with loading is higher for a laminate with a short delamination of 10 mm compared to a laminate with long delamination of 26 mm. This implies that a laminate with a short delamination is more sensitive to loading as a small change in load can significantly influence the delamination growth rate which is governed by changes to G. However, it is important to note that $\epsilon_{cr}$ is significantly higher for short delamination at 0.0050, compared to 0.0007 for long delamination. Consequently, the laminate with short delamination stores more membrane energy due to compression before G starts to increase, which may contribute to a steeper slope in the increase of G. In the next subsection, insights into growth characteristics of delamination under fatigue and quasi-static loading conditions have been discussed.

### 3.2.2. Delamination growth characteristics

The Fig. 7 presents two graphs in which Fig. 7(a) presents the insights into growth characteristics under fatigue loading. This has been done by plotting the normalised energy release rate, $G/G_c$ against delamination length ($l$), for various end-shortening strains, $\epsilon_0$. Fig. 7(b) presents the corresponding decrease in axial stiffness with $l$. Initially Fig. 7(a) has been discussed which is followed later by Fig. 7(b). To determine the conditions for stable and unstable growth under fatigue, it is sufficient to examine the trend of $G/G_c$ under the maximum strain associated with fatigue loading. Therefore as also discussed previously, in Fig. 7, $\epsilon_0$ represents maximum end-shortening strain and G represents $G_{max}$ for fatigue related discussions.

As already demonstrated in subsection 3.2.1, G begins to increase once the thin sublamine undergoes buckling. Similarly, in this case, when $\epsilon_0$ is held constant and $l$ is increased, the thin sublamine buckles at a critical delamination length, at which point G starts to increase above zero. As expected, for higher $\epsilon_0$, there is a reduction in critical length at which G begins to increase. The fatigue growth behaviour depends on the trend of G following the buckling of the thin sublamine. Based on this, the design space can be categorised into two growth regimes: fatigue and static, as illustrated in Fig. 7(a). In the fatigue growth regime, fatigue growth occurs as long as $G > G_{th}$ and is governed by Paris law, which predicts the growth rate as a function of G. Maximum strain allowable for which the growth is expected to remain in this regime is for $\epsilon_a = 0.00217$ (black curve in Fig. 7(a)). In this regime, if G increases with delamination length $l$, the growth is considered unstable; conversely, if G decreases with $l$, the growth is considered stable (Chai et al., 1981; Köllner et al., 2020). To understand this better, consider a laminate with an initial delamination length of 20 mm subjected to $\epsilon_0 = 0.00217$. Initially, $G/G_c$ is 0.841, reaching its peak value of 1.0 at a delamination length of 27 mm (points marked in Fig. 7(a)). Beyond this point, $G/G_c$ steadily decreases, reaching a minimum value of 0.598 at a delamination length of 98 mm. As G initially increases with delamination length from 20 mm to 27 mm, the growth rate is expected to accelerate, indicating unstable growth. Beyond this point, G gradually decreases with increasing delamination length from 27 mm to 98 mm, indicating a stable growth. A further discussion on delamination growth rates is presented in the following subsection on imperfection sensitivity. However, in the fatigue growth regime, the delamination is not expected to grow abruptly to a larger length, as is the case in the static growth regime. In the static growth regime, the growth is governed by condition for quasi-static loading, which requires $G \ge G_c$ for growth to occur. As an example, for $\epsilon_0 = 0.00250$, the growth is expected to be abrupt between delamination lengths of 18 mm and 58 mm, where $G/G_c \ge 1.0$, but beyond 58 mm, the growth is expected to transition back into the fatigue growth regime. For $\epsilon_0 = 0.00285$, once the growth enters static regime at delamination length of 16 mm,

<figure><img src="image_12.png" /><figcaption>**Fig. 7.** For thick laminate: (a) Normalised energy release rate ($G/G_c$) with respect to delamination length ($l$) at different end-shortening strains ($\epsilon_0$). Insights into growth characteristics under fatigue loading are presented. (b) Normalised stiffness degradation $((F/\Delta)_{\text{Normalised}})$ with respect to $l$.</figcaption></figure>

it does not transition back to the fatigue regime as shown in Fig. 7(a). The similar trend is observed for $\epsilon_0$ values of greater than 0.00285.

Overall, it can be concluded that for the growth to remain in the fatigue regime, $\epsilon_0$ should not exceed $\epsilon_a = 0.00217$ for the selected laminate configurations and loading conditions. To slow down the growth rates further, $\epsilon_0$ of less than $\epsilon_a$ can be selected as the allowable limit. Another important factor to consider as the delamination length increases is the degradation of the axial stiffness of the laminate. The Fig. 7(b) illustrates the normalised axial stiffness as a function of delamination length. The axial stiffness is determined by dividing the axial reaction force (F) in intact part (part (3)) of the delaminated laminate by the applied constant end-shortening displacement $\Delta$, and it is normalised relative to the stiffness of a laminate without delamination. The plots are presented for two end-shortening strains from within the fatigue growth regime. The results indicate that while stiffness decreases with increasing delamination length, the reduction is not very significant. This is attributed to the thick sublamine, which has a much higher stiffness compared to the thin sublamine, acting as the primary load-bearing element. Due to negligible bending of the thick sublamine, the stiffness degradation remains minimal, with only about a 5% reduction even at a delamination length of 95 mm. All the insights presented above are crucial for selecting the appropriate parameters for fatigue analyses planned in future studies and, consequently, for applying the concept of 'slow-growth' damage tolerance to the problem of delamination buckling.

Next, the growth characteristics under quasi-static loading conditions in Fig. 8 are presented and related to fatigue growth characteristics in Fig. 7(a). In Fig. 8, the end-shortening strain ($\epsilon_0$) is plotted on x-axis versus the delamination length ($l$) on the y-axis. For a given $l$, as $\epsilon_0$ increases, it first reaches the buckling strain of the thin sublamine (part (1) of laminate), $\epsilon_{cr}$, indicated by the dashed blue line in Fig. 8. As expected, $\epsilon_{cr}$ decreases as the $l$ increases. With further increase beyond $\epsilon_{cr}$, the quasi-static growth criterion ($G = G_c$) is met at the failure strain $\epsilon_f$, represented by the solid black line. Based on the delamination
<figure><img src="image_13.png" /><figcaption>**Fig. 8.** For thick laminate: End-shortening strain ($\epsilon_0$) on x-axis) at the buckling ($\epsilon_{cr}$) of the thin sublaminate and failure ($\epsilon_f$) (when the delamination growth criterion is met) as a function of delamination length ($l$) (on y-axis) under quasi-static loading conditions. Insights into growth characteristics under quasi-static loading are presented.</figcaption></figure>

length $l$, the $\epsilon_f$ curve is divided into three parts, each representing a distinct growth regime. These regimes: unstable, unstable → stable and stable are depicted in different shades of grey. Unstable growth occurs when $l \le l_u$ (Fig. 8), where, upon reaching $\epsilon_f$ for a given $l$, the delamination grows abruptly to nearly the full length of the laminate. For instance, for the $l = 14$ mm, the $\epsilon_f = 0.00320$ which is greater than $\epsilon_f$ of all the delamination lengths from 15 mm and 90 mm. In the unstable→ stable regime, growth initially occurs in an unstable manner but transitions later to stable manner, requiring an increase in $\epsilon_0$ for further progression. This regime occurs when $l_u < l < l_s$ (Fig. 8), such as for $l = 22$ mm. In this case, $\epsilon_f = 0.00223$, which is greater than the $\epsilon_f$ for all delamination lengths from 23 mm to 35 mm but beyond $l = 35$ mm, $\epsilon_0$ has to be increased for the growth to progress. The final stable regime occurs when $l \ge l_s$ (Fig. 8) where $l_s = 26$ mm. For instance for $l = 26$ mm, $\epsilon_f = 0.00217$ and subsequently $\epsilon_f$ increases with $l$, which implies that $\epsilon_0$ has to be increased for the growth to progress, which in turn means that this regime exhibits stable growth. A key point to be noted is that in both Fig. 7(a) and Fig. 8, the abrupt growth occurs when $G > G_c$. This is the reason, the minimum value of $\epsilon_{fmin} = 0.00217$ in Fig. 8 corresponds to the maximum allowable strain ($\epsilon_a$) in the fatigue regime depicted in Fig. 7(a). With this insight, a few quasi-static FE simulations can be conducted for relevant delamination lengths of 14 to 35 mm to determine the allowable strain ($\epsilon_a$) for growth to remain in fatigue regime. The findings from the quasi-static FE analysis later will offer deeper insights into the advantages of the adopted analytical method and highlight potential future directions for this research. More details on this will follow towards the end of the paper in subsection 3.4. It is known that geometric imperfection can influence the buckling behaviour which subsequently also influence the growth behaviour. Therefore, the next subsection presents a discussion on the influence of imperfection on the growth characteristics.

### 3.2.3. Imperfection sensitivity

Three magnitudes of imperfections have been studied which are defined by the maximum perturbation ($q_0$) at the centre of part (1) of the laminate. The values of $q_0$, which have been analysed are 0.001, 0.01 and 0.02. Fig. 9 illustrates the influence of imperfections on normalised energy release rate $G/G_c$, which is plotted against the delamination length ($l$) for two different values of end-shortening strain $\epsilon_0$, 0.00189 and 0.00217. The previous subsection demonstrated that $\epsilon_0 = 0.00217$ corresponds to maximum allowable strain ($\epsilon_a$), for growth to stay within the fatigue regime. However, this applies only when $q_0 = 0.01$ or lower. If however, $q_0$ is higher – for instance, 0.01 or 0.02, then $G$ values may slightly enter the static growth regime, as shown in Fig. 9, which can have a minor influence on the strain allowables.

Similarly, imperfection magnitudes may also impact delamination growth rates under fatigue, consequently affecting the number of cycles

<figure><img src="image_14.png" /><figcaption>**Fig. 9.** Normalised energy release rate ($G/G_c$) versus the delamination length ($l$). Influence of imperfection magnitude ($q_0$) has been presented.</figcaption></figure>

required for the delamination to grow by a prescribed length. This has been analysed for $\epsilon_0 = 0.00189$ which lies within the fatigue regime. Fig. 9 shows that for this $\epsilon_0$ the $G$ values are also marginally affected due to different values of $q_0$. For the analysis of growth rate, it is assumed that the laminate initially has a 20 mm delamination, and fatigue growth is evaluated until the delamination length reaches 60 mm. The Mode I Paris law, derived from a double cantilever beam (DCB) test, is used in this analysis, which is based on a study conducted by Dávila et al. (2020). The Paris Law used in this analysis has the form $\frac{dl}{dN} = 820.0 G_{max}^{7/4}$. As discussed in previous subsections the $G$ shown in Fig. 9 corresponds to maximum strain in context of fatigue loading and thus represents $G_{max}$. The parameters in Paris law above correspond to ratio of minimum to maximum energy release rate of 0.1. The results of this analysis are presented in Fig. 10. In the previous graphs (Figs. 7 to 9), the normalised energy release rate ($G/G_c$) and failure strain ($\epsilon_f$) have been presented for different delamination lengths ($l$) without considering delamination growth. However, in Fig. 10, where delamination growth is considered, the delamination length is denoted as $l_g$.

Fig. 10(a) shows the delamination growth rate ($\Delta l_g/\Delta N$) plotted against delamination length ($l_g$), while Fig. 10(b) displays $l_g$ as a function of the number of cycles ($N$). The $\Delta l_g/\Delta N$ follows the same trend as $G$ which is shown in Fig. 9. For all the imperfection magnitudes, $G$ value initially increases, as delamination length increases from 20 mm to 30 mm. After reaching the peak at 30 mm delamination, $G$ declines gradually as the length increases to 60 mm. Similarly, $\Delta l_g/\Delta N$ initially increases with the growth of $l$ from 20 to 30 mm and subsequently decreases as $l$ grows from 30 to 60 mm. The main observation in the results is that imperfection significantly influences the $\Delta l_g/\Delta N$, which in turn influences the number of cycles required for the growth of delamination (Fig. 10(b)). The number of cycles required for the delamination to grow from 20 to 60 mm is around 114000 for $q_0 = 0.001$. For $q_0 = 0.01$ the number of cycles decreases by 16.7% reducing to 95000. For $q_0 = 0.02$, the number of cycles decrease by 32.4%, resulting in 77000 cycles. This study therefore demonstrates that imperfections may significantly influence the postbuckling growth behaviour of delaminations under fatigue loading, and should therefore be considered in the design processes. However, in all other sections of this study, an imperfection magnitude of $q_0 = 0.001$ has been adopted. It is important to note that selecting higher imperfection magnitudes will not affect the qualitative outcomes of the results.

### 3.3. Discussions on thin laminate

The thin laminate refers to a laminate with the layup of $[0_3/ -45_3/45_3/90_3]_{2S}$. The results discussed in first two subsections are for this thin laminate with near-surface through-width delamination, located between the final $0_3$ layer and the penultimate $-45_3$ layer. However, the final subsection presents an in-depth discussion on how the through-thickness location of delamination influences the buckling
<figure><img src="image_15.png" /><figcaption>**Fig. 10.** Fatigue delamination growth evaluation at $\epsilon_0 = 0.00189$ for different imperfection magnitudes ($q_0$): (a) Delamination growth rate ($\Delta l_g/\Delta N$) versus delamination length ($l_g$). (b) Delamination length ($l_r$) versus the number of cycles ($N$).</figcaption></figure>

response and growth behaviour. In all the subsections, the maximum imperfection magnitude at the centre of part (1), $q_0$, is defined as 0.001. Unlike in thick laminates, the imperfection sensitivity study is not presented for thin laminates, as the results are qualitatively similar to those observed in the thick laminate.

### 3.3.1. Insight into sensitivity of delamination growth to load

Similar to the thick laminate, initially insights into sensitivity of delamination growth to fatigue loading has been presented by comparing a short delamination length of 10 mm with a long one of 26 mm. This has been done again by examining the rate of increase of $G$ with applied end shortening strain ($\epsilon_0$) under quasi-static loading conditions. The results for delamination length of 10 mm are presented in Fig. 11, while those for delamination length of 26 mm are shown in Fig. 12. Similar to the case of the thick laminate, Fig. 11(a) and Fig. 12(a) feature two y-axes, both plotted against $\epsilon_0$ on the x-axis. The membrane energy ($U_{Mem}^{(1)}$) and bending energy ($U_{Bend}^{(1)}$) contributions in part (1) of the laminate are shown on the primary y-axis, while $G/G_c$ is plotted on the secondary y-axis. In Fig. 11(b) and Fig. 12(b) out-of-plane deformations at the centre of part (1) ($w^{(1)}$) and part (2) ($w^{(2)}$) of the laminates have been illustrated.

The results exhibit many similarities to the thick laminate. Initially, as $\epsilon_0$ increases in both laminates with delamination lengths of 10 mm and 26 mm, only $U_{Mem}^{(1)}$ increases, reflecting axial shortening of the laminates. This behaviour continues until the buckling strain ($\epsilon_{cr}$) of the thin sublamine, or part (1) of the laminate, is reached. Beyond this point, $w^{(1)}$ begins to increase, leading to a rise in $U_{Bend}^{(1)}$, while $U_{Mem}^{(1)}$ starts to plateau. After buckling, $G/G_c$ begins to increase for both the delamination lengths. Similar to the thick laminate, $G/G_c$ exhibits a much steeper increase with loading for the 10 mm delamination length compared to the 26 mm length. The $\epsilon_f$ values of 0.00535 for a delamination length of 10 mm and 0.00193 for a delamination

<figure><img src="image_16.png" /><figcaption>**Fig. 11.** For thin laminate and delamination length of 10 mm: (a) Membrane ($U_{Mem}^{(1)}$) and bending ($U_{Bend}^{(1)}$) energies in the thin sublamine and normalised energy release rate ($G/G_c$) with respect to end shortening strain ($\epsilon_0$). (b) Out-of-plane displacements of thin ($w^{(1)}$) and thick ($w^{(2)}$) sublaminates with respect to $\epsilon_0$.</figcaption></figure>

length of 26 mm are 1.5% and 11.05% lower than the corresponding values for the thick laminate. The difference is due to slightly higher values of $w^{(2)}$, with -0.005 mm for the 10 mm delamination length and -0.046 mm for the 26 mm length. This occurs because part (2), or the thick sublamine, is significantly thinner in this case compared to the thick laminate discussed in Subsection 3.2.1. Here, $\epsilon_f$ is compared for only two delamination lengths, however, the next subsection provides a more comprehensive overview across various delamination lengths.

The insights into growth sensitivity concerning fatigue loading remain consistent with those discussed earlier for thick laminates. The results indicate that for short delaminations, the steeper increase in $G$ with $\epsilon_0$ makes it more sensitive to loading compared to long delaminations. The next subsection discusses the growth characteristics of delamination under fatigue and quasi-static loading conditions, comparing these results with those previously reported for thick laminates.

### 3.3.2. Delamination growth characteristics

Fig. 13 presents two graphs with Fig. 13(a) depicting the growth characteristics under fatigue loading and Fig. 13(b) showing the degradation of normalised axial stiffness ($F/\Delta$) Normalised of the laminate with respect to delamination length ($l$). Similar graphs for thick laminates have been presented in Section 3.2.2, which can be referred for basic description on them. Similar to the thick laminate, Fig. 13(a) illustrates two growth regimes: the fatigue and the static one. It should be noted again that in Fig. 13, $\epsilon_0$ represents maximum end-shortening strain and $G$ represents $G_{max}$ for fatigue related discussions. The fatigue growth behaviour is influenced by the trend of $G$ following the buckling of the thin sublamine. As mentioned in Subsection 3.2.2, in the fatigue regime, the growth occurs when $G > G_{th}$ and follows Paris law, which describes the growth rate as a function of $G$. Compared to the thick laminate, thin laminates exhibit higher peak values of $G$, followed
<figure><img src="image_17.png" /><figcaption>**Fig. 12.** For thin laminate and delamination length of 26 mm: (a) Membrane ($U_{Mem}^{(1)}$) and bending ($U_{Bend}^{(1)}$) energies in the thin sublaminate and normalised energy release rate ($G/G_c$) with respect to end shortening strain ($\epsilon_0$). (b) Out-of-plane displacements of thin ($w^{(1)}$) and thick ($w^{(2)}$) sublaminates with respect to $\epsilon_0$.</figcaption></figure>

by a more pronounced decline after the peak for different values of $\epsilon_0$. For example, at $\epsilon_0 = 0.00160$, the peak value of $G/G_c$ is 0.700 at delamination length of 30 mm, which subsequently decreases with delamination length to 0.239 at 98 mm delamination. On the other hand for the thick laminate, at the same $\epsilon_0 = 0.00160$, the peak value is 0.541 for delamination length of 30 mm, which then decreases to 0.340 for 98 mm length (refer to Fig. 7(a)).

Since the peaks of G are higher for thin laminate, maximum strain allowable for which growth is expected to remain in fatigue regime ($\epsilon_a$) also decreases compared to thick laminate. For the thin laminate, the value of $\epsilon_a$ is 0.00192 which is 11.5% lower than the value for the thick laminate at 0.00217. Furthermore, since the post peak decline of G is more pronounced for thin laminates, delamination growth after the peak is expected to slow down more under fatigue compared to thick laminates. In the static growth regime, since G exceeds $G_c$, the growth is governed by condition for quasi-static loading, where growth occurs when $G \ge G_c$. Similar trend of higher peaks and more pronounced declines thereafter is observed in static regime. For example, for $\epsilon_0 = 0.00285$, in case of thin laminate, once G enters static regime, it transitions back to fatigue regime due to steep decline later. On the other hand, in case of thick laminate (Fig. 7(a)), for same $\epsilon_0$, G stays in the static regime after entering it.

Overall, it can be summarised that maximum allowable strain ($\epsilon_a$), for the growth to remain in fatigue regime decreases with the reduction in laminate thickness. However, the fatigue growth rate is expected to slow down more after the peak values of G is reached. Another factor to consider is stiffness degradation with delamination growth. Since thin laminates exhibit higher out-of-plane deformation of thick sublaminate, they are expected to undergo higher stiffness degradation with delamination length. As mentioned above, the normalised axial stiffness degradation for the thin laminate is depicted in Fig. 13(b)

<figure><img src="image_18.png" /><figcaption>**Fig. 13.** For thin laminate: (a) Normalised energy release rate ($G/G_c$) with respect to delamination length ($l$) at different end-shortening strains ($\epsilon_0$). Insights into growth characteristics under fatigue loading are presented. (b) Normalised stiffness degradation $((F/\Delta)_{Normalised})$ with respect to $l$.</figcaption></figure>

as a function of delamination length. It is similarly presented for two $\epsilon_0$ values within the fatigue growth regime. The results show that the thin sublaminate exhibits significantly higher stiffness degradation than the thick laminate, as expected. At a delamination length of 95 mm, the stiffness degradation is approximately 15%, which is three times higher than that of the thick laminate. Therefore, stiffness degradation must be considered in the design process to evaluate the load-carrying capacity of the laminate as delamination progresses. All these insights are essential for future fatigue studies and for applying the principles of 'slow-growth' damage tolerance to address delamination buckling challenges.

Fig. 14 illustrates the growth characteristics under quasi-static loading conditions. A similar graph has been discussed for the thick laminate in subsection 3.2.2, which can be referred to for the basic description. Three distinct growth regimes are depicted: unstable, unstable→ stable and then stable. Unstable growth occurs to nearly full length of the laminate for $l \le l_u$, as the $\epsilon_f$ corresponding to this value $l$ exceeds failure strains of all the subsequent delamination lengths. Similar to thick laminate, the unstable→ stable regime occurs for $l_u < l < l_s$. In this regime, the growth is initially expected to be unstable but later transitions to stable, which implies that $\epsilon_0$ has to be increased for growth to progress. As an example, consider a laminate with delamination of 22 mm, for which $\epsilon_f = 0.00202$. This value exceeds all the $\epsilon_f$ values until delamination of 40 mm, beyond which $\epsilon_0$ has to be increased for the growth to progress. In thin laminate unstable→ stable regime is broader compared to thick laminate shown in Fig. 7(b). This is because, this regime is bounded between the minimum value of $\epsilon_{fmin}$ and the value of $\epsilon_f$ corresponding to final delamination length of 90 mm. For thin laminate $\epsilon_{fmin}$ is 0.00192 which is lower than $\epsilon_{fmin}$ for thick laminate at 0.00217. However, for delamination length of 90 mm, the value of $\epsilon_f$ is higher for the thin laminate at 0.00336, compared to 0.00277 for the thick laminate. This reversal in the trend, where $\epsilon_f$ is higher for the thin laminate compared to the thick laminate, is observed starting from a delamination length of 56 mm and beyond. This trend reversal is attributed to the steeper decline of G post the peak, as
<figure><img src="image_19.png" /><figcaption>**Fig. 14.** For thin laminate: End-shortening strain ($\epsilon_0$) on x-axis at the buckling ($\epsilon_{cr}$) of the thin sublaminate and failure ($\epsilon_f$) when the delamination growth criterion is met) as a function of delamination length ($l$) on y-axis under quasi-static loading conditions. Insights into growth characteristics under quasi-static loading are presented.</figcaption></figure>

shown in Fig. 13(a). Consequently, the applied $\epsilon_0$ must be increased for larger delamination lengths for $G$ to reach $G_c$. The final stable regime occurs for $l \ge l_s$, where $\epsilon_0$ has to be increased for growth to progress. Another similarity to thick laminate is that, the minimum value of $\epsilon_{fmin} = 0.00192$ in Fig. 14 corresponds to the maximum allowable strain ($\epsilon_a$) in the fatigue growth regime depicted in Fig. 13(a). It can thus be summarised that, for thin laminates, the unstable → stable regime is broader compared to thick laminates, and failure strains are higher for larger delamination lengths starting from 56 mm.

Since it is unclear to what extent the results from the specific configurations presented so far can be generalised to arbitrary through-thickness delamination locations, the following subsection addresses this aspect in detail for the case of thin laminate. Finally, based on the insights presented here and in subsection 3.2.2, a quasi-static FE study is proposed to aid in the verification of the reported outcomes of the present work. A comprehensive discussion of the results from the quasi-static FE study is presented in subsection 3.4.

### 3.3.3. Influence of delamination depth

The influence of through-thickness location of the delamination has been presented for the case of thin laminate. The configurations examined so far in this study featured thick-to-thin sublaminate thickness ratios of 47:1 for the thick laminate and 15:1 for the thin laminate. In the configurations explored in this section, this ratio has been further reduced to 14:2, 13:3, and 11:5. These studies have shown that the laminate can also exhibit a global, or closed, buckling mode response, which may lead to penetration between the thin and thick sublaminates at the delaminated interface. This has been handled by applying a constraint condition ($w^{(1)}(x_1 = l/2) - w^{(2)}(x_2 = l/2) \ge 0.0$, where the term in brackets represents the difference in out-of-plane displacement between the thin and thick sublaminates at the centre of the laminate. This condition is employed following the approach of Köllner (2021), by augmenting the total strain energy ($U$) with the term $\lambda(w^{(1)}(x_1 = l/2) - w^{(2)}(x_2 = l/2))$, where $\lambda$ is the Lagrange multiplier representing the contact force. Applying the variational principle of the Total Potential Energy ($\delta\Pi = 0$), then results in an additional equation $\frac{\partial\Pi}{\partial\lambda} = 0.0$ to determine the unknown $\lambda$.

The original configuration of thin laminate (Configuration 1), which is already presented in subsections 3.3.1 and 3.3.2 can be represented as: $[[0_3/ - 45_3/45_3/90_3]_2/90_3/45_3/ - 45_3/0_3/90_3/45_3/ - 45_3//0_3]$ where ‘//’ denotes the through-thickness location of delamination. In addition, the following configurations have been investigated:

* Configuration 2: $[[0_3/ - 45_3/45_3/90_3]_2/90_3/45_3/ - 45_3/0_3/90_3/45_3// - 45_3/0_3]$, where the delamination is located between the layers 14 and 15.

* Configuration 3: $[[0_3/ - 45_3/45_3/90_3]_2/90_3/45_3/ - 45_3/0_3/90_3//45_3/ - 45_3/0_3]$, where the delamination is located between the layers 13 and 14.

* Configuration 4: $[[0_3/ - 45_3/45_3/90_3]_2/90_3/45_3/ - 45_3//0_3/90_3/45_3/ - 45_3/0_3]$, where the delamination is located between the layers 11 and 12.

Similar to configuration 1 of thin laminate, the analyses have been conducted starting with delamination length ($l$) of 10 mm for the additional configurations mentioned above. Initially the results of configuration 2 have been presented. The Fig. 15(a) presents the buckling modes observed in this configuration for delamination lengths of 10 mm, 16 mm and 20 mm. This is done by plotting the out-of-plane displacements at the centre of thin ($w^{(1)}$) and thick ($w^{(2)}$) sublaminates along the x-axis, against the applied end-shortening strain ($\epsilon_0$) on the y-axis. The Fig. 15(b) shows the corresponding increase in energy release rate $G$ on the y-axis with respect to the $\epsilon_0$ on the x-axis. For the 10 mm delamination case, the laminate exhibits a closed global buckling mode response from the outset. This behaviour contrasts with that of the original thin laminate configuration with a 10 mm delamination, as shown in Fig. 11(b). In that case, only an opening mode response is observed after the limit point strain of 0.005 is reached. Since the buckling response in the present configuration is in the closed mode, delamination failure is likely governed by Mode II fracture. Therefore, the analysis has been carried out until $\epsilon_0 = 0.0115$ when $G_{IIC}$ is reached as shown in Fig. 15(b). At this point, the out-of-plane displacements at the centre of both thin and thick sublaminates are 0.04 mm. Similar behaviour is observed up to a 12 mm delamination; however, the laminate response changes beyond this length. For example, in the case of a 16 mm delamination, the laminate initially exhibits a closed global buckling mode. This response persists up to the limit point $\epsilon_0 = 0.0078$, beyond which an opening mode buckling response is observed and the two sublaminates exhibit out-of-plane displacements in opposite directions (see Fig. 15(a)). In this case, the limit point corresponds to a value of $G$ that lies between $G_{IIC}$ and $G_{IIC}$ and $G$ shoots up due to change in buckling mode (see Fig. 15(b)). For 20 mm delamination, similar buckling response is observed, however $G$ at limit point is lower than $G_{IIC}$.

For longer delaminations, the overall buckling response remains similar. However, the limit point at which the mode transition occurs decreases as the delamination length increases. For example, Fig. 16(a) below illustrates the buckling responses for delamination lengths of 30 mm, 40 mm, and 74 mm. A closer view is also provided in Fig. 16(b) to enhance the visualisation of the displacement $w^{(2)}$. The Fig. 16 also indicates the points at which mode I failure occurs for different delamination lengths following open mode buckling.

Similar to configuration 1 of thin laminate, Fig. 17 depicts the growth characteristics under fatigue loading by plotting $G/G_c$ versus the delamination length $l$. In this figure $G_c$ corresponds to mode I failure. It has been shown that maximum strain allowable $\epsilon_a$ for which the growth is expected to be in the fatigue regime is 0.00182 which has reduced compared to 0.00192 for original configuration (Fig. 13(a)). But the peak of $G/G_c$ corresponding to $\epsilon_a$ occurs at a higher delamination length of 48 mm, compared to 28 mm in original configuration (Fig. 13(a)).

Next, configuration 3 is discussed briefly. In this case, the buckling responses are similar to those observed in configuration 2. Therefore, for the sake of brevity, the buckling response is not shown for this configuration. In this case, for delamination length ranging from 10 mm to 22 mm, a closed global buckling mode is observed up to strains corresponding to mode II failure ($G_{IIC}$). In contrast, the configuration 2 exhibited this behaviour only up to a delamination length of 12 mm. For delamination lengths beyond 22 mm, the buckling response is observed to transition from closed to opening mode upon reaching a limit point. This behaviour is consistent with that seen in the case of Configuration 2. At a delamination length of 34 mm, the value
<figure><img src="image_20.png" /><figcaption>**Fig. 15.** For configuration 2: (a) Buckling response for laminate with short delaminations illustrated by out-of-plane displacements of thin ($w^{(1)}$) and thick ($w^{(2)}$) sublaminates with respect to $\epsilon_0$. (b) Corresponding increase in energy release rate $G$ with respect to $\epsilon_0$.</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>**Fig. 16.** For configuration 2: (a) Buckling response for laminate with long delaminations illustrated by out-of-plane displacements of thin ($w^{(1)}$) and thick ($w^{(2)}$) sublaminates with respect to $\epsilon_0$. (b) Closer view for better visualisation of $w^{(2)}$.</figcaption></figure>

of $G$ at the limit point falls below $G_{IC}$. This has been observed to occur at the delamination length of 20 mm for configuration 2. Beyond

<figure><img src="image_22.png" /><figcaption>**Fig. 17.** For configuration 2: Normalised energy release rate ($G/G_c$) with respect to delamination length ($l$) at different end-shortening strains ($\epsilon_0$). Insights into growth characteristics under fatigue loading are presented.</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>**Fig. 18.** For configuration 3: Normalised energy release rate ($G/G_c$) with respect to delamination length ($l$) at different end-shortening strains ($\epsilon_0$). Insights into growth characteristics under fatigue loading are presented.</figcaption></figure>

this point, the limit point strain continues to decrease with increasing delamination length, following a similar trend. The fatigue growth characteristics plotted in Fig. 18 below shows that $\epsilon_a$ for which the growth is expected to be in fatigue regime increases to 0.00194. This value exceeds both the 0.00182 observed for configuration 2 and the 0.00192 for the configuration 1. The peak value corresponding to $\epsilon_a$ is observed at a much higher delamination length of 76 mm.

Finally, the results of configuration 4 have been presented. For this configuration, the buckling responses corresponding to various delamination lengths have been illustrated in Fig. 19(a) below. In Fig. 19(b) a closer view is shown for better visualisation of the $w^{(2)}$. Like the previous configurations, for short delaminations closed global mode response is observed until $\epsilon_0$ corresponding to mode II failure ($G_{IIc}$) is reached. For configuration 4 this is observed for delamination lengths up to 40 mm which is much higher than corresponding values for previous configurations. For longer delaminations, such as 56 mm, an initial closed global mode response is observed, which transitions to an opening mode after reaching the limit point at $\epsilon_0 = 0.0067$ but there is a difference observed in this case compared to previous configurations. For configuration 4, after reaching the limit point, the thick sublamine moves in the direction of thin sublamine, but $w^{(1)}$ is significantly higher than $w^{(2)}$. Although the sublaminates move in the same direction, the buckling response remains in the open mode due to differences in their out-of-plane displacement. Also, for this configuration, 56 mm is the least delamination length for which $G$ at limit point is lower than $G_{IC}$. For higher delamination lengths identical buckling response has been observed as shown in Fig. 19(a) and (b) for the lengths of 66 mm and 80 mm. The figure also shows points corresponding to mode II failure for delamination length of 40 mm and mode I failure corresponding to delamination lengths of 66 mm and 80 mm. For the 56 mm delamination length, the mode I failure point is not shown, as the value of $G$ abruptly increases to significantly higher value after mode transition. The fatigue growth characteristics plotted in Fig. 20 below shows that $\epsilon_a$ for which the growth is expected to be
<figure><img src="image_24.png" /><figcaption>**Fig. 19.** For configuration 4: (a) Buckling response for laminate with delaminations illustrated by out-of-plane displacements of thin ($w^{(1)}$) and thick ($w^{(2)}$) sublaminates with respect to $\epsilon_0$. (b) Closer view for better visualisation of $w^{(2)}$.</figcaption></figure>

<figure><img src="image_25.png" /><figcaption>**Fig. 20.** For configuration 4: Normalised energy release rate ($G/G_c$) with respect to delamination length ($l$) at different end-shortening strains ($\epsilon_0$). Insights into growth characteristics under fatigue loading are presented.</figcaption></figure>

in fatigue regime increases further to 0.0025, but the value of 1.0 of $G/G_c$ is observed almost at the end of laminate at 98 mm.

In summary, as the through thickness location of delamination increases, the laminate initially exhibits closed global mode deformation, which persists up to mode II failure in the case of short delaminations. For longer delaminations, although the initial response is also closed mode, the laminate transitions to an open mode after reaching a limit point. This shift in buckling response from closed to open mode has also been reported by Köllner et al. (2020) for a cross-ply laminate. Additionally, it is observed that for shallower delaminations, the sublaminates move in opposite directions after the limit point, whereas for deeper delaminations, they move in the same direction. Finally, it is shown that for delamination growth to remain in fatigue regime, the threshold $\epsilon_a$ initially decreases with increasing delamination depth, but it later increases as the depth increases further.

### 3.4. Quasi-static failure strain curves from FE model

This subsection aims to verify the key outcomes of the analytical model presented in Subsections 3.2.2 and 3.3.2 through nonlinear FE

analyses. The key difference from the verification example in Section 2.2 is that this section considers two distinct failure criteria: one is a more conservative mode I (opening) criterion, similar to that in Section 2.2, while the other is a mode dependent criterion. The FE analysis has been conducted in ABAQUS, with the plies modelled using solid elements (C3D8I) and the interfaces represented by cohesive elements. Delamination growth has been simulated using a fracture mode mixity-dependent cohesive zone model (CZM) based on Turon et al. (2018). The CZM parameters for the unidirectional graphite/epoxy composite (IM7/8552) have been adopted from Joosten et al. (2018) and are presented in Table 3 below. In the more conservative case of the mode I-based failure criterion, the strengths and critical energy release rates for mode II (shearing) and mode III (tearing) are set equal to those of mode I. However, it should be noted that for graphite/epoxy materials like the one used in this study, the fracture process zone (FPZ) is sufficiently small, making linear elastic fracture mechanics (LEFM) applicable. This suggests that delamination growth under variable fracture mode loading is governed by energy release rate, with interlaminar strength having no influence on the energy dissipated. The FE model is described in detail in Venkat et al. (2025).

Fig. 21(a) presents the FE predictions, showing the variation of $\epsilon_f$ with $l$, where the selected thick and thin laminates are compared for the two discussed failure criteria. In Fig. 21(a), the minimum failure strain, $\epsilon_{fmin}$, corresponding to the four curves is also highlighted (points a. to d.). Fig. 21(b) illustrates the fracture mode mixity at failure as a function of $l$, comparing the selected thick and thin laminates. This plot is shown only for calculations based on the mode dependent failure criterion, as the fracture mode influences the $\epsilon_f$ in these cases. Fracture mode mixity is defined as $\frac{G_{II}+G_{III}}{G_I+G_{II}+G_{III}}$, where $G_I$, $G_{II}$ and $G_{III}$ represent the energy release rates associated with mode I, mode II and mode III delamination growth, respectively. These $G$ terms are determined by calculating the respective areas under the traction-separation in CZM law in FE analyses. Although mode III is included in this formulation, its contribution to the problem of a delaminated strut under compression is negligible compared to the other modes. A fracture mode mixity of 0.0 represents pure mode I, while a value of 1.0 represents pure mode II. Values between 0.0 and 1.0 correspond to a mixed mode fracture. In Fig. 21, the marker points represent the delamination lengths that have been analysed and these are connected by dotted lines to visualise the expected trend in between them.

Initially results in Fig. 21(a) based on mode I criterion have been discussed and compared to the analytical results for both the thick ($\epsilon_f$ in Fig. 8) and thin laminates ($\epsilon_f$ in Fig. 14). It should be noted that in Figs. 8 and 14, the x-axis represents $\epsilon_0$ since the graphs depict both $\epsilon_{cr}$ and $\epsilon_f$. In contrast, Fig. 21(a) displays only $\epsilon_f$. Consistent with analytical results, the FE results show that as the delamination length increases from 12 mm to 26 mm, $\epsilon_f$ decreases from 0.00348 to minimum value ($\epsilon_{fmin}$) of 0.00205 for the thick laminate and from 0.00336 to $\epsilon_{fmin}$ of 0.00192 for thin laminate. The $\epsilon_{fmin}$ values show excellent agreement with the analytical results, which predict values of 0.00217 for thick laminates and 0.00192 for thin laminates. Beyond 26 mm, the values for both thick and thin laminates continue to increase, with the values for the thick laminate remaining higher up to 50 mm. Subsequently, the $\epsilon_f$ values of the thin laminate exceed those of the thick laminate and at a delamination length of 90 mm, the values are 0.0031 for the thick laminate and 0.0036 for the thin laminate. These findings are again in very good agreement with the analytical results, where the $\epsilon_f$ values at 90 mm delamination length have been found to be 0.0028 for the thick laminate and 0.0034 for the thin laminate. Overall, the results from the analytical model closely align with the FE predictions based on the mode I failure criterion, making it a useful tool for gaining initial insights into the problem.

In Fig. 21(a) the FE predictions based on the mode dependent failure criterion similarly exhibits several qualitative similarities to the analytical results. Initially, as the delamination length increases from 12 mm to 18 mm, $\epsilon_f$ decreases from 0.0036 to minimum value ($\epsilon_{fmin}$)
<table><caption>**Table 3**
Interface cohesive zone properties of IM7/8552 unidirectional graphite/epoxy composite from Joosten et al. (2018).</caption><thead><tr><th colspan="5">Cohesive zone properties</th></tr><tr><th>τ<sub>3</sub><sup>o</sup> [MPa]</th><th>τ<sub>2</sub><sup>o</sup> = τ<sub>2</sub><sup>o</sup> = τ<sub>sh</sub><sup>o</sup> [MPa]</th><th>G<sub>lc</sub> [N/mm]</th><th>G<sub>IIc</sub> = G<sub>IIIc</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>

<figure><img src="image_26.png" /><figcaption>**Fig. 21.** Results of FE analyses: (a) Failure strain curve ($\epsilon_f$) (on x-axis) for different delamination lengths ($l$) (on y-axis). The selected thick and thin laminates are compared for different failure criteria. (b) Fracture mode mixity (on x-axis) at failure for different $l$. The selected thick and thin laminates are compared for mode dependent failure criterion.</figcaption></figure>

of 0.0027 for the thick laminate and from 0.0035 to $\epsilon_{fmin}$ of 0.00248 for thin laminate. It is to be noted that, as with analytical results, the $\epsilon_{fmin}$ of thick laminate is higher than that of thin laminate. Subsequently, as the delamination length is increased from 18 mm to 50 mm, $\epsilon_f$ continues to increase for both the laminate types and up to 50 mm the $\epsilon_f$ values for the thick laminate remain higher than those for the thin laminate. After this point, the trend reverses, with $\epsilon_f$ values for the thin laminate exceeding those of the thick laminate. At the delamination length of 90 mm, the $\epsilon_f$ for thin laminate is 0.0068, which is higher than the corresponding value of thick laminate of 0.00552. These findings confirm that these FE predictions demonstrate clear qualitative alignment with the analytical results. Nonetheless, noticeable quantitative discrepancies between the analytical and these predictions based on mode dependent failure criterion are observed from the present study. One of the main differences is related to the $\epsilon_{fmin}$ values, which are under-predicted by the analytical approach. This is because the fracture criterion is based on mode I $G_c$, a conservative assumption, as the presented analytical method does not account for the fracture mode mixity at this stage. However, it is clear from Fig. 21(b) that at the delamination length of 18 mm, corresponding to $\epsilon_{fmin}$, the fracture

mode mixities are 0.55 for the thick laminate and 0.50 for the thin laminate, which represent mixed mode fracture growth. Fig. 21(b) shows that, in general, for smaller delamination lengths, growth occurs under mixed-mode conditions, with the mode II component increasing as the delamination length increases. For delamination lengths beyond 50 mm, the growth takes place mainly due to mode II as fracture mode mixity is close to 1.0. It is due to the increasing fracture mode mixity that $\epsilon_f$ becomes significantly higher compared to analytical results as the delamination length increases from 18 mm to 90 mm. These findings are in line with an earlier study by Kardomateas et al. (1995), in which the increase in fracture mode mixity with increasing delamination length have been reported. Considering the same, this approach of including the influence of fracture mode mixity into the analytical formulations could be considered as a future scope of this research.

Overall, it can be concluded that the adopted analytical model provided valuable insights into the problem of delamination buckling, yielding results that qualitatively match the FE findings. The analytical model allows for the determination of the delamination length range at which $\epsilon_{fmin}$ is likely to occur. As discussed in the previous section, this $\epsilon_{fmin}$ value represents the maximum allowable strain that should not be exceeded to ensure growth remains within the fatigue regime. Such insights are essential for significantly reducing the number of FE analyses that need to be performed. On the other hand, the results from the analytical model are conservative, as it does not account for the fracture mode mixity and uses mode I $G_c$ as the fracture criterion, potentially paving the way for its further refinements in future.

## 4. Conclusions and outlook

Delaminations in composite structures can severely compromise their mechanical performance by increasing the risk of premature buckling and promoting further growth. In this study, a Rayleigh-Ritz based analytical model proposed by Köllner et al. (2020) has been used to gain insights into the growth behaviour of delaminations for compressive fatigue loading. The growth behaviour has been analysed and compared for examples of thick and thin laminates with generic quasi-isotropic (QI) layups and a near-surface through-width delamination. The key findings of the analytical study are summarised below.

* Insights into the sensitivity of delamination growth to fatigue loading have been presented by comparing laminates with initial delamination lengths of 10 mm and 26 mm. For both thick and thin laminates, the analysis demonstrated that the increase in $G$ with $\epsilon_0$ is steeper for short delamination (10 mm) compared to long (26 mm) one. This indicates that short delaminations are more sensitive to loading.

* The growth characteristics under fatigue loading have been examined for both thick and thin laminates, identifying two distinct growth regimes: fatigue and static. The maximum allowable strain limit ($\epsilon_a$) required to maintain growth within the fatigue regime has been defined and it has been shown that $\epsilon_a$ for thin laminate is 11.5% lower than for thick laminate. Furthermore, the trend of $G$ which governs the growth rate, demonstrated differences between thick and thin laminate. The growth characteristics under quasi-static loading has also been investigated and linked with fatigue growth characteristics. This has been done by showing that minimum failure strain ($\epsilon_{fmain}$) under quasi-static loading corresponds to $\epsilon_a$.
* The analysis of geometric imperfection sensitivity has been conducted by evaluating the growth rate using Paris Law. The results revealed that the growth rate follows the trend of G and is significantly influenced by imperfections. Consequently, the number of cycles required for the delamination to grow by a certain length is also affected due to imperfections.

* The influence of delamination depth on the results is also discussed. It is shown that the buckling response varies with delamination depth, and that it also affects the maximum allowable strain limit, $\epsilon_a$, needed to ensure that growth remains within the fatigue regime. It has been shown that $\epsilon_a$ initially decreases with delamination depth, but subsequently increases as the depth continues to increase.

Finally, the outcomes of the analytical model have been verified by determining the $\epsilon_f$ as a function of delamination length ($l$) using geometrically nonlinear FE analysis. It has been shown that FE predictions based on mode I failure criterion have very good agreement with outcomes based on analytical model. For mode dependent failure criterion, although the FE predictions qualitatively match the analytical results, a quantitative discrepancy is noticeable. This can be attributed to the influence of fracture mode mixity on the FE results, which is not accounted for in the analytical results presented in this study. The influence of the fracture mode mixity can, however, be incorporated in future studies using the formulations presented in Köllner (2021), which examine buckling-driven delamination growth in plates with a circular embedded delamination. This formulation accounts for the fracture mode mixity through an approach based on the work of Davidson et al. (2000), where crack-tip forces and moments are used to calculate both the total energy release rate and the individual mode-dependent contributions. These contributions can then be used to determine the fracture mode mixity. In conclusion, this study, utilising the energy based analytical approach, provides valuable insights into the growth characteristics in the delamination buckling problem under fatigue loading. These insights can aid in determining critical parameters for future FE based fatigue studies, such as $\epsilon_a$, which must not be exceeded to prevent abrupt delamination growth and to ensure that growth stays within the fatigue regime. Such insights will consequently enable designers to apply the less conservative 'slow-growth' damage tolerance concept to address delamination buckling.

### CRediT authorship contribution statement

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

### Declaration of AI-assisted technologies in the writing process

During the preparation of this work, the authors used AI tools to improve the language and grammar checks. After using these tools, the authors reviewed and edited the content as needed and take full responsibility for the content of the publication.

### 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. Basic formulations of the analytical framework

The additional derivations in this section follow from Köllner et al. (2020). The laminate is clamped on both ends. These end boundary conditions along with continuity conditions at the interface between various laminate parts are given below

$$ 
\begin{aligned}
\omega^{(3)}(x_3 = 0) &= 0, & \omega^{(3)}(x_3 = sl) &= \omega^{(k)}(x_k = 0), \\
\omega^{(3)}(x_3 = L) &= 0, & \omega^{(3)}(x_3 = sl + l) &= \omega^{(k)}(x_k = l), \\
(\omega^{(3)})'(x_3 = 0) &= 0, & (\omega^{(3)})'(x_3 = sl) &= (\omega^{(2)})'(x_k = 0), \\
(\omega^{(3)})'(x_3 = L) &= 0, & (\omega^{(3)})'(x_3 = sl + l) &= (\omega^{(k)})'(x_k = l)
\end{aligned} 
\quad (A.1) $$

where $k = 1$ to 2 represents the first two parts of the laminate while terms corresponding to part (3) are specified in the equation. The postbuckling displacements of the three parts of the laminate are approximated using Rayleigh–Ritz method by enforcing the boundary and continuity conditions. The displacements for parts 1 and 2 are as follows

$$ \omega^{(k)} = q_k \sin^2\left(\frac{\pi x_k}{l}\right) + C_0 x_k^2 + C_1 x_k^2 + C_2 x_k + C_3, \quad \text{where } k = 1, 2 \text{ and coefficients} $$

$$ 
\begin{aligned}
C_0 &= \frac{q_3}{2l^3} \left( \frac{L}{\pi} \tan\left(\frac{\pi sl}{L}\right) - \frac{L}{\pi \cos\left(\frac{\pi sl}{L}\right)} \sin\left(\frac{\pi (sl + l)}{L}\right) \right) \\
C_1 &= \frac{1}{2} \left( -\frac{2q_3}{L} - 3C_0 l \right), \quad C_2 = q_3, \quad C_3 = \frac{q_3 L}{2\pi} \tan\left(\frac{\pi sl}{L}\right)
\end{aligned} 
\quad (A.2) $$

while for part 3

$$ \omega^{(3)} = q_3 \frac{L}{2\pi \cos\left(\frac{\pi sl}{L}\right) \sin\left(\frac{\pi sl}{L}\right)} \sin^2\left(\frac{\pi x_3}{L}\right) \quad (A.3) $$

The axial shortening in the three parts of laminate, based on nonlinear kinematic approach in Hunt et al. (2004) are shown below

$$ 
\begin{aligned}
u^{(1)} &= q_4 - \frac{1}{2} \int_0^l \left( \frac{\partial w^{(1)}}{\partial x_1} \right)^2 dx_1 - 2 \left( \frac{1-a}{2} \right) t q_3, \\
u^{(2)} &= q_4 - \frac{1}{2} \int_0^l \left( \frac{\partial w^{(2)}}{\partial x_2} \right)^2 dx_2 + 2 \left( \frac{a}{2} \right) t q_3,
\end{aligned} 
\quad (A.4) $$

Finally, with axial shortening and postbuckling displacement, the in-plane strains and curvatures are defined in each laminate part as follows

$$ \epsilon_0^{(k)} = \frac{u^{(k)}}{l}, \quad \epsilon_0^{(3)} = \frac{u^{(3)}}{sl}, \quad \kappa_0^{(k)} = \frac{\partial^2 w^{(k)}}{\partial x_k^2}, \quad \kappa_0^{(3)} = \frac{\partial^2 w^{(3)}}{\partial x_3^2} \quad (A.5) $$

where $k = 1$ to 2.

### Data availability

The necessary data required to replicate the findings are available in the manuscript. If any additional data are required, they may be made available upon request.
References

Ameri, B., Moradi, M., Mohammadi, B., Salimi-Majd, D., 2020. Investigation of nonlinear post-buckling delamination in curved laminated composite panels via cohesive zone model. Thin-Walled Struct. 154, 106797. http://dx.doi.org/10.1016/j.tws.2020.106797.

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

Bisagni, C., Brambilla, P., Davila, C.G., 2013. Modeling delamination in postbuckled composite structures under static and fatigue loads. In: Proceedings of International SAMPE Technical Conference. Long Beach, California, United States, pp. 1035–1049.

Brunner, A.J., 2008. Experimental study of delamination in cross-ply laminates. In: Delamination Behaviour of Composites. In: Woodhead Publishing Series in Composites Science and Engineering, Elsevier, pp. 281–309. http://dx.doi.org/10.1533/9781845694821.3.281.

Chai, H., Babcock, C.D., 1985. Two-dimensional modelling of compressive failure in delaminated laminates. J. Compos. Mater. 19 (1), 67–98. http://dx.doi.org/10.1177/002199838501900105.

Chai, H., Babcock, C.D., Knauss, W.G., 1981. One dimensional modelling of failure in laminated plates by delamination buckling. Int. J. Solids Struct. 17 (11), 1069–1083. http://dx.doi.org/10.1016/0020-7683(81)90014-7.

Davidson, Barry D., Yu, LiJie, Hu, Hurang, 2000. Determination of energy release rate and mode mix in three-dimensional layered structures using plate theory. Int. J. Fract. 105, 81–105. http://dx.doi.org/10.1023/a:1007672131026.

Dávila, C.G., Rose, C.A., Murri, G.B., Jackson, W.C., Johnston, W.M., 2020. Evaluation of fatigue damage accumulation functions for delamination initiation and propagation. Technical report, Langley Research Center.

Harper, P.W., Hallett, S.R., 2008. Cohesive zone length in numerical simulations of composite delamination. Eng. Fract. Mech. 75 (16), 4774–4792. http://dx.doi.org/10.1016/j.engfracmech.2008.06.004.

Huang, H., Kardomateas, G.A., 1998. Buckling of orthotropic beam-plates with multiple central delaminations. Int. J. Solids Struct. 35 (13), 1355–1362. http://dx.doi.org/10.1016/s0020-7683(97)00110-8.

Hunt, G.W., Hu, B., Butler, R., Almond, D.P., Wright, J.E., 2004. Nonlinear modeling of delaminated struts. AIAA J. 42 (11), 2364–2372. http://dx.doi.org/10.2514/1.5981.

Joosten, M.W., Yang, Q.D., Blacklock, M., Wang, C.H., 2018. A cohesive network approach for modelling fibre and matrix damage in composite laminates. Compos. Struct. 206, 658–667. http://dx.doi.org/10.1016/j.compstruct.2018.08.087.

Kardomateas, G., Pelegri, A.A., Malik, B., 1995. Growth of internal delaminations under cyclic compression in composite plates. J. Mech. Phys. Solids 43 (6), 847–866. http://dx.doi.org/10.1016/0022-5096(95)00012-8.

Kharazi, M., Ovesy, H.R., 2008. Postbuckling behavior of composite plates with through-the-width delaminations. Thin-Walled Struct. 46 (7), 939–946. http://dx.doi.org/10.1016/j.tws.2008.01.005.

Kharazi, M., Ovesy, H.R., Taghizadeh, M., 2010. Buckling of the composite laminates containing through-the-width delaminations using different plate theories. Compos. Struct. 92 (5), 1176–1183. http://dx.doi.org/10.1016/j.compstruct.2009.10.019.

Köllner, Anton, 2021. Predicting buckling-driven delamination propagation in composite laminates: An analytical modelling approach. Compos. Struct. 266, 113776. http://dx.doi.org/10.1016/j.compstruct.2021.113776.

Köllner, A., Jungnickel, R., Völlmecke, C., 2016. Delamination growth in buckled composite struts. Int. J. Fract. 202 (2), 261–269. http://dx.doi.org/10.1007/s10704-016-0158-y.

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

Köllner, A., Völlmecke, C., 2018. Post-buckling behaviour and delamination growth characteristics of delaminated composite plates. Compos. Struct. 203, 777–788. http://dx.doi.org/10.1016/j.compstruct.2018.03.010.

Liu, P.F., Hou, S.J., Chu, J.K., Hu, X.Y., Zhou, C.L., Liu, Y.L., Zheng, J.Y., Zhao, A., Yan, L., 2011. Finite element analysis of postbuckling and delamination of composite laminates using virtual crack closure technique. Compos. Struct. 93 (6), 1549–1560. http://dx.doi.org/10.1016/j.compstruct.2010.12.006.

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

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

Mohammadi, B., Shahabi, F., 2016. 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. 152, 88–104. http://dx.doi.org/10.1016/j.engfracmech.2015.04.005.

Nikrad, S.F., Keypoursangari, S., Asadi, H., Akbarzadeh, A.H., Chen, Z.T., 2016. Computational study on compressive instability of composite plates with off-center delaminations. Comput. Methods Appl. Mech. Engrg. 310, 429–459. http://dx.doi.org/10.1016/j.cma.2016.07.021.

Ovesy, H.R., Kharazi, M., 2011. Compressional stability behavior of composite plates with through-the-width and embedded delaminations by using first order shear deformation theory. Comput. Struct. 89 (19), 1829–1839. http://dx.doi.org/10.1016/j.compstruc.2010.10.016, Civil-Comp.

Pascoe, J.-A., 2020. Slow-growth damage tolerance for fatigue after impact in FRP composites: Why current research won't get us there. Procedia Struct. Integr. 28, 726–733.

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

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

Sheinman, I., Kardomateas, G.A., Pelegri, A.A., 1998. Delamination growth during pre- and post-buckling phases of delaminated composite laminates. Int. J. Solids Struct. 35 (1-2), 19–31. http://dx.doi.org/10.1016/s0020-7683(97)00058-9.

Sheinman, I., Soffer, M., 1991. Post-buckling analysis of composite delaminated beams. Int. J. Solids Struct. 27 (5), 639–646. http://dx.doi.org/10.1016/0020-7683(91)90218-5.

Simitzes, G.J., Sallam, S., Yin, W.L., 1985. Effect of delamination of axially loaded homogeneous laminated plates. AIAA J. 23 (9), 1437–1444. http://dx.doi.org/10.2514/3.9104.

Sridharan, S., 2008. Delamination Behaviour of Composites. Woodhead Publishing, pp. 1–762.

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

Venkat, S.S., Scheffler, S., Anilkumar, P.M., Baranger, E., Rolfes, R., 2025. Effect of delamination defects on buckling and growth characteristics in composite laminates. Compos. Struct. 119245. http://dx.doi.org/10.1016/j.compstruct.2025.119245.

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

Zhang, G., Hu, Y., Yan, B., Tong, M., Wang, F., 2024. Buckling and post-buckling analysis of composite stiffened panels: A ten-year review (2014–2023). Thin-Walled Struct. 112525. http://dx.doi.org/10.1016/j.tws.2024.112525.