## Abstract

An efficient analytical model is described which predicts the value of compressive strain below which buckle-driven propagation of delaminations in aerospace composites will not occur. An extension of this efficient strip model which accounts for propagation transverse to the direction of applied compression is derived. In order to provide validation for the strip model a number of laminates were artificially delaminated producing a range of thin anisotropic sub-laminates made up of 0°, ±45° and 90° plies that displayed varied buckling and delamination propagation phenomena. These laminates were subsequently subject to experimental compression testing and nonlinear finite element analysis (FEA) using cohesive elements. Comparison of strip model results with those from experiments indicates that the model can conservatively predict the strain at which propagation occurs to within 10 per cent of experimental values provided (i) the thin-film assumption made in the modelling methodology holds and (ii) full elastic coupling effects do not play a significant role in the post-buckling of the sub-laminate. With such provision, the model was more accurate and produced fewer non-conservative results than FEA. The accuracy and efficiency of the model make it well suited to application in optimum ply-stacking algorithms to maximize laminate strength.

## 1. Introduction

Carbon fibre-reinforced plastics (CFRPs) have the potential to radically reduce the weight of any structure in which they are employed. This is a consequence of their favourable strength and stiffness properties and the way in which these can be tailored to give different directional properties. However, even in aircraft such as the Boeing 787 and the Airbus A350 this weight saving may not be fully exploited. There are a number of reasons for this; one of the most significant being barely visible impact damage (BVID), a topic which has been the subject of substantial research effort in the field of composites for the last 30 years. When layered structures are impacted their structural properties can be severely reduced. Often this is a consequence of the formation of hidden delaminations. This problem becomes acute when damaged structures are placed under compressive loads which tend to cause thin, outer sub-laminates to buckle locally, creating a mechanism for delamination propagation and premature failure. Hence aerospace regulations specify that such laminates must tolerate BVID without failure for ultimate levels of load, where ultimate is usually 50 per cent above limit load. This typically reduces allowable strains to around 4500 μstrain. In some applications, this gives component strengths that are less than half of the equivalent values for aluminium.

The analysis of delamination propagation following buckling deformation of a sub-laminate is complicated by the occurrence of a mix of peeling (mode I), shearing (mode II) and tearing (mode III) actions along the delamination perimeter. Hence, nonlinear finite element analyses (FEAs) employing interface elements or virtual crack closure techniques are the methods usually adopted (e.g. [1–4]). However, such methods are complex, requiring large computational effort, and do not lend themselves to understanding the mechanisms that promote or reduce damage tolerance across a range of laminates. It is therefore crucial that simplified methods are developed to capture such mechanisms and thus provide engineering tools to improve damage tolerance during optimum laminate design.

The authors have previously derived a simple, semi-numerical, fracture mechanics model to predict critical values of applied strain at which such buckle-driven propagation of delamination occurs [5]. This work extended the principles of strain energy release rate (SERR) considered previously [6–8] to the case of a two-dimensional, anisotropic plate containing a delamination. In this paper, the authors present a new extension to this extremely efficient strip model [5]. The model assumes simplified components of bending and membrane (strain) energy in the post-buckled, thin sub-laminate created by the delamination and is based on the concept of establishing an equivalent mode I approximation of the actual, mixed-mode SERR. The current paper clarifies the underlying methodology of the previous work [5] and extends it by improving the accuracy of the strip model for transverse growth, perpendicular to the loading direction. This is accomplished by considering the direction of maximum axial stiffness in the thin sub-laminate in either the loading or transverse directions.

The new strip model methodology and its simplifying assumptions are validated against the results of a series of new experiments and FEAs on laminates containing a single artificial circular delamination producing different types of anisotropy in the thin, buckled sub-laminate. The experiments were monitored using a high-speed digital image correlation (DIC) system which produces plots of full-field deformation and principal strains so that the occurrence of critical buckling as well as the onset, mode and stability of propagation can be pin-pointed.

## 2. Modelling methods

### (a) Strip model: the concept of equivalent energy release

#### (i) Initial post-buckling

Hunt *et al.* [9] showed that a thin sub-laminate created by near-surface delamination of an isotropic strut placed under compressive load can buckle to open the delamination. A transition from delamination opening to closing occurs for deeper delaminations where the point of transition defines the critical depth. This depth is dependent on the effective length of the strut compared with the length of the delamination. When the strut has an effective length that is twice the delamination length, the critical depth occurs at about 25 per cent of total thickness.

For the laminated plates that follow, and for most delaminations that have been caused by impact damage, the delaminations from which propagation will occur are at depths less than 25 per cent of total thickness. Furthermore, the anti-buckling guides used in the tests described below create effective lengths less than twice the delamination length. Hence these delaminations open under compressive load, and the failure is triggered by this opening (sub-laminate buckling).

#### (ii) One-dimensional propagation

The proposed model assumes that all post-buckled energy in the thin sub-laminate available for propagation of the delamination is in the form of bending energy, i.e. there is no membrane energy available after buckling, meaning that there is no release of energy associated with post-buckled stiffness.

Initially, a one-dimensional, thin-film sub-laminate buckling above a flat substrate (figure 1), is developed so that the principles can later be extended to the more complicated two-dimensional case. Thin-film conditions are assumed so that no bending energy is released from the substrate during propagation.

Applying a sinusoidal buckling displacement
2.1
The bending energy *U*_{B} in the buckled sub-laminate, which has bending rigidity *EI*, is given by
2.2
The end-shortening of the sub-laminate is taken as
2.3
Hence the potential energy *V* of the force *P* applied to the thin sub-laminate is
2.4
For equilibrium, according to the principle of minimum potential energy,
2.5
Substituting equations (2.2)–(2.4) into equation (2.5), we obtain the critical buckling load of the sub-laminate
2.6
For unit width and thickness *t* of the sub-laminate, the critical strain is then
2.7
If we ignore any post-buckling stiffness and assume that the end-shortening of equation (2.3) is equal to the applied end-shortening following buckling, then
2.8
From equation (2.2), the bending energy in the buckle can be re-written as
2.9

The rate at which bending energy is released as the delamination propagates is
2.10
Assuming that there is no post-buckled stiffness in the sub-laminate, we add the membrane energy released by the small unbuckled element *δL* in figure 1 as the delamination propagates,
2.11
Similarly, the membrane energy retained by the sub-laminate is given by
2.12
The difference of equations (2.11) and (2.12) gives the energy available for release and therefore the rate of membrane energy released during propagation of element *δL*,
2.13
This allows for the fact that only the membrane energy above the critical strain is available for release from the un-buckled element.

Adding together the bending and membrane energies released as a result of propagation of the buckled sub-laminate, we obtain the following SERR:
2.14
Note that the final terms of equations (2.10) and (2.13) cancel, meaning that the result is unchanged if we assume no variation of *ε*^{C} with *L* [5]. Equation (2.14) was obtained by Chai *et al*. [6] and requires the assumption that there is no post-buckled stiffness in the thin sub-laminate. If a further assumption is made that any mode II contribution is ignored and a pure mode I peeling action is assumed at the tip of the delamination, a mode I propagation criterion can be applied. The mode II contribution may actually be significant but, since the energy required for mode I fracture is always significantly less than that required for mode II fracture, the mode I criterion is conservative. Hence, propagation is assumed to occur when *G*_{IC}, the critical SERR for the material in which the delamination occurs, is reached (*G*=*G*_{IC}), at which point the applied strain reaches its threshold propagation value (*ε*=*ε*_{th}), and
2.15
Re-arranging for the threshold strain gives
2.16

#### (iii) Two-dimensional buckling and propagation

The two-dimensional sub-laminate is assumed to buckle as a thin film above a flat substrate (figure 2). In the case considered here, we assume that out-of-plane rotation and out-of-plane displacement are prevented along the circular perimeter of the delamination and that the sub-laminate is subject to forces arising from compatibility with the full laminate (substrate). These forces are calculated from the full laminate stiffness matrix in the following way. Assuming zero curvature in a balanced and symmetric parent laminate (displaying neither extension/shear coupling nor in-plane/out-of-plane coupling), the loads acting on the thin sub-laminate {*N*}_{SL} are determined by obtaining the strain {*ε*}_{L} of the full parent laminate when unit axial strain is applied. {*N*}_{SL} is then calculated, by assuming compatibility of strain, from
2.17
i.e.
2.18
where [*A*]_{SL} is the in-plane (membrane) stiffness matrix of the thin sub-laminate and *ν* is the effective Poisson ratio of the full laminate.

Note that, although uni-axial strain is applied to the full laminate (shear strain, *γ*_{xy}=0), equation (2.17) may result in bi-axial load and shear being applied to the thin sub-laminate. Note also that equation (2.17) assumes that positive strains and loads are compressive.

In order to calculate the critical value of the applied axial strain *ε*^{C} an extremely efficient, infinite strip software, Viconopt [10], is used in the results that follow. The values of critical strain could also be determined using finite element software, as is shown later for comparison, but the use of the infinite strip model is much more efficient. This approach divides the circular sub-laminate into a number of strips of equal width while applying compatibility and equilibrium along the connecting nodal lines. The differential equations representing periodic buckling deformation of these strips are then solved exactly to develop a transcendental eigenvalue problem [11]. Points at which these nodal lines intersect with the circular delamination boundary are restraining points where the sinusoidal wavelengths are combined to match the boundary using energy minimization principles [12]. For the results presented later, six equal width strips are used with 12 constrained nodes at the junction of these strips and the circular boundary (see [13]). The thin, two-dimensional sub-laminate can be unbalanced and asymmetric, which can give rise to fully populated matrices for in-plane membrane stiffness [*A*]_{SL}, out-of-plane bending stiffness [*D*]_{SL} and coupling stiffness [*B*]_{SL}. However, Viconopt buckling analysis is fully general, and can analyse such laminates.

The post-buckling behaviour of the thin sub-laminate is now schematically described in order to explain the actual mechanism by which strain energy is released as the delamination propagates in two dimensions. This will become useful when we make simplifying assumptions about the SERR. Considering the two orthogonal strips AA and BB shown in figure 2, we see that the post-buckled compressive strip (figure 2*b*) is restrained by downward forces representing the effect of orthogonal tensile strips (e.g. figure 2*d*), that are being stretched to accommodate the post-buckling deformation. The combined effects give rise to mixed-mode conditions all along the delamination front [14] with transverse propagation being dominated by mode II effects. As above in the one-dimensional case, we isolate the mode I (peeling) effect by assuming that all strain energy available for propagation is in the form of (i) bending energy in the thin sub-laminate and (ii) membrane energy (above *ε*^{C}) in the unbuckled perimeter of the delamination. This energy is then released by a peeling mechanism, so that the (transverse) orthogonal strips are not assumed to develop any tensile components, and hence an *equivalent* mechanism for mode I propagation is produced. Removing these tensile components means that the sub-laminate has no membrane stiffness following buckling, i.e. radial compression is assumed to act, as illustrated in figure 2*c*,*e*. The basis for these assumptions is the fact that bending energy is released as if the sub-laminate is a strut, i.e. in proportion to the inverse of the delaminated length squared. In actual fact, the mode I contribution must be less than this since the two-dimensional shape has post-buckled stiffness owing to the development of tensile stresses in the lateral direction. Hence, less bending energy than is represented by equation (2.14) is actually released. However, it is expedient, with regard to producing a simple model, that ignoring these nonlinear two-dimensional effects appears to have little influence on the prediction of threshold propagation strain for thin sub-laminates, as shown by comparison with the experimental and FEA results that follow. This appears to be because the rate of energy release as propagation occurs is independent of the absolute value of strain energy in the post-buckled state.

Hence, an equivalent mode I energy criterion is introduced, replacing the need for a mixed-mode expression and nonlinear post-buckling analysis. The strain *ε*_{th} at which propagation occurs is as above, i.e.
2.19
where *A*_{nn,SL} is the axial stiffness of the sub-laminate and
or
In the former case (*n*=1), the longitudinal stiffness of the sub-laminate dominates and so propagation is assumed to initiate in the *x*-direction and SERR is calculated for this direction. In the latter case (*n*=2), propagation initiates in the *y*-direction, we assume that pre-buckling compressive strain is applied in the longitudinal *x*-direction in order to find *ε*^{C} and that post-buckled compressive strain is applied in the transverse *y*-direction in order to develop the equivalent mode I energy release in that direction owing to radial compression. This assumes that the sub-laminate buckles as a result of the actual applied strain but all subsequent strain energy is developed in the direction of greater axial stiffness.

### (b) Finite element analysis using cohesive elements

In order to provide validation for the strip model and to give further insight into the supporting experimental validation, the nonlinear finite element program Abaqus [15] was applied to predict the strains required for sub-laminate buckling and delamination growth in the artificially delaminated laminates considered in the next section. The laminate was separated into three regions, with the interface zone, which is populated with cohesive elements, being sandwiched between the thin sub-laminate and the thick base laminate (figure 3). Four-noded shell elements (S4) were used in the sub-laminate and base to account for any through-thickness shear deformation that may arise in the test coupons. In addition, the area around the delamination front was finely meshed to ensure solution convergence and to accurately capture delamination initiation and growth. Note that the three layers had the same mesh design and density (see figure 4), and that the nodes on the bottom edge had clamped boundary conditions. Compression was applied as a uniform displacement in the *x*-direction along the top edge of the mesh where all other displacements and rotations were fixed. Moreover, the nodes between the circular edge of the anti-buckling guide described below and the four straight edges were constrained to prevent out-of-plane displacement, as is the function of the anti-buckling guide in the physical experiments described in §3*a*.

In order to prevent overlap of the sub-laminate and base, thus ensuring feasible solutions, a surface-to-surface contact constraint was applied in the delamination zone. Eight-noded three-dimensional cohesive elements were used at the interface between the base and sub-laminates, in order to predict delamination propagation. The cohesive layer was composed of zero thickness volumetric elements, which have the following definitions.

Stiffness is defined to relate the stress to the relative displacement at the interface. This parameter was chosen such that the value was large enough to ensure that the sub-laminate and base were well bonded prior to damage growth and small enough not to cause convergence problems. Three stiffness parameters can be defined in terms of peeling, shearing and tearing components, i.e. modes I, II and III. The critical interlaminar material strength for each fracture mode was defined using a traction-separation law. For pure mode I loading, material damage occurs once the interfacial normal stress reached its maximum interlaminar tensile strength, following which the stiffness was gradually reduced to zero. Under pure mode I, II or III loading, the onset of damage at the interface can be determined simply by comparing the stress components with their respective maximum. However, under mixed-mode loading, damage onset may occur before any of the individual stress components reaches their respective maximum value. Therefore, a damage initiation criterion was defined to predict this mixed-mode fracture. In this study, it was assumed that delamination initiation could be predicted using the quadratic failure criterion [1]. A damage propagation criterion, based on fracture energy, was defined to predict delamination growth under mixed-mode conditions. Here the energy-based Benzeggagh–Kenane [16] law, which is the simplest to implement and most commonly used, was employed to predict delamination growth. The total critical energy release rate *G*_{TC} is given by
2.20
where *η* is an empirical parameter derived from mixed-mode tests, *G*_{IC}, *G*_{IIC} and *G*_{IIIC} are critical energy release rates for pure modes I, II and III, respectively, and *G*_{I}, *G*_{II} and *G*_{III} are the corresponding values determined from analysis. According to the criterion, damage growth occurs when the total energy *G*_{I}+*G*_{II}+*G*_{III}≥*G*_{TC}. Within the Abaqus model, there are two stages involved in propagation. The first stage is delamination initiation at which the nodes at the delamination front start to separate; the second stage is delamination propagation in which entire cohesive elements fail and are removed.

## 3. Results

In order to test the applicability of the strip model, a series of experimental and FEA test cases were considered. Experimental validation consisted of a series of compression tests on laminates constructed from unidirectional CFRP layers. Each laminate was designed to display local sub-laminate buckling and propagation behaviour in order to provide benchmark results for the strip model. This was achieved by employing non-stick circles of polytetrafluoroethylene (PTFE) at a single interface within each laminate. A range of examples (see table 1) were created by varying the degree of anisotropy displayed by the thin sub-laminate above the PTFE insert, which is dependent on the number and orientation of plies making up the sub-laminate.

Note that positive fibre orientations are defined as being clockwise with respect to the applied load direction (*x*-axis). Similarly, for the FEA test cases, a model was constructed using shell and cohesive elements, as explained above. This FEA model incorporated an anti-buckling guide in order to mimic the physical experiments as closely as possible. Eight so-called quasi-isotropic laminates, each having equal percentages of CFRP plies with fibres aligned at 0°, +45°, −45° or 90° to the direction of applied load, were manufactured from 0.25 mm thick Hexcel T700GC/M21 pre-preg CFRP layers, which have material properties *E*_{11}=136 GPa, *E*_{22}=8.9 GPa, *G*_{12}=4.5 GPa, *ν*_{12}=0.35 and *G*_{1C}=550 J m^{−2}. Sequences of layers that make up the laminates are given in table 1. Note that the full laminate stacking sequences are balanced and symmetric and thus negate the occurrence of the coupling effects described in the text above equation (2.17). During manufacture, PTFE circles 0.0125 mm thick and 39 mm in diameter were placed between selected layers in the laminates to produce areas where adhesion would be prevented, i.e. artificial delaminations. Through-thickness positions of these inserts are given in table 1 by indicating the thin sub-laminate created. The test coupons will subsequently be referred to using the sub-laminate stacking sequence. Note that coupons with sub-laminates containing more than one 90° ply have artificial delaminations at deeper interfaces than the other coupons as modelling indicated that a [90_{2}] sub-laminate would have a delamination buckling strain that was too high to induce failure as a consequence of propagation following local sub-laminate buckling.

### (a) Experimental test set-up

Prior to compression testing, laminates were cut into coupons 210 mm long by 100 mm wide by 4 mm thick and placed in a compression fixture with an integrated circular anti-buckling guide of internal diameter 85 mm in order to prevent buckling outside this circle (see figure 5). Tests consisted of applying axial compression under displacement control at 0.1 mm min^{−1} until local propagation and/or global failure occurred. The faces of the coupons nearest to the PTFE inserts were covered in a random speckle pattern to allow buckling modes and failure sequences to be visualized (following post-processing) using a DIC system. This system employs a pair of stereo cameras to produce plots of out-of-plane displacement or variation in surface strain in the principal load direction relative to a reference image taken under zero load.

As verification of ‘far-field’ strain for the DIC system and to ensure that specimens were correctly aligned, strains were recorded throughout the tests by two pairs of vertically aligned back-to-back strain gauges. For each experiment, one pair of gauges was located inside but as close as possible to the anti-buckling guide on a vertical line through the centre of the circular PTFE insert, and similarly the second pair was located inside but as close as possible to the anti-buckling guide on a horizontal line through the centre of the PTFE insert. See figure 6 for a schematic of strain gauge placement for the [45_{2}] laminate.

### (b) Experimental results

Results from the experimental tests are given in table 2. Experimental strains for local buckling and propagation are determined from a combination of load–strain plots (figure 6) and DIC images (figure 5). Unless noted otherwise, for all DIC images colours indicate changing out-of-plane displacement from an unloaded reference state. The DIC images and load–strain plots were linked by the load output of the compression test machine. The occurrence of propagation was detected by fitting a circle of diameter equal to the PTFE insert (using the internal anti-buckling guide diameter as a reference) around the fully formed buckle on a DIC image then visually determining when the local buckle had spread outside of this circle (figure 5*b*,*c*). Using the LED load readout, propagation was correlated with the average of the four strain gauge readings on the respective load–strain plot to determine an experimental strain for initiation of propagation. The average of the four strain gauge readings was used as it accounts for losses in stiffness during the test.

#### (i) Local buckling

Local buckling manifests itself as a change in gradient of the strain gauge curves (e.g. figure 6) and as an area with a high number of contours on the DIC images, implying a rapidly changing out-of-plane displacement (figure 5). However, these changes are often subtle, as can be seen from figure 5*a*; hence, for a number of results in table 2 a range of buckling strains is given. The lower value coincides with when a buckle may have occurred on the DIC images and the upper value with when a buckle has definitely formed (contrast figure 5*a*,*b*). In general, local buckling mode shapes for sub-laminates with no 0° ply had a smaller wavelength in the load direction, resulting in a laterally elongated local buckle that did not fill the artificially delaminated area (figure 7*b*). Note that the random pattern outside the local buckle seen in figure 5*a* is partially a result of the system being unable to resolve the very small deflection changes occurring at low levels of load in this area and partly a result of the non-smooth surface of the sample which is a consequence of the surface finish.

During some experiments (indicated by superscript ‘a’ in table 2), despite the PTFE insert, some adhesion remained between the sub- and base laminates. Hence, the sub-laminates for these coupons were not released until sufficient load had been applied to overcome this adhesion, resulting in an increased buckling strain. For the [45_{2}], [0/90_{2}], [0_{2}/90] and [0/45] sub-laminates, sufficient load was applied to release the sub-laminate prior to propagation, after which loading was removed and the test reset. Thus buckling occurred for these tests without adhesion. In the case of the [0/90_{2}] coupon, no adhesion between the PTFE insert and sub-laminate occurred. During initial loading of this coupon (table 2), a loud noise was mistaken for a release of the sub-laminate and the test was reset. DIC images taken during the initial test identified that the loud noise was due to a sudden change in local buckling mode shape (figure 7*a*,*b*). This so-called mode jump caused intra-ply cracking (between fibres in a ply), which resulted in local buckling in the second test bypassing the first buckling mode shape and instead exhibiting the second mode shape at a lower level of strain. This accounts for the discrepancy between buckling strains in table 2. Changes in local buckling mode shape were also noted in the [0_{2}/90] and [0/45] sub-laminates, although in both these cases local buckling mode jumps coincided with a global buckling/bending event. The occurrence of global bending or buckling phenomena can be seen as wider contours in figure 5 in comparison with the thin contours and hence steep gradient changes that represent local buckling. Figure 8*a*,*b* shows that as the load increases significant global deformation occurred for the [0_{2}/90] sub-laminate.

#### (ii) Propagation

Propagation was deemed to have occurred once the buckle spread outside a 39 mm diameter circle centred at the point of greatest out-of-plane deflection (contrast figure 5*b*,*c*). This propagation was correlated with the corresponding load–strain plot, where it was generally seen as a discontinuity in the strain gauge curves (see figure 6), in which multiple propagation events are evident. However, as for local buckling, if growth occurred slowly it was not always possible to pinpoint when propagation occurred. Slow growth can occur as a stiffness change rather than a discontinuity on the load–strain plots. Thus, for some cases in table 1 a range of propagation strains is given. The direction of propagation was influenced to a greater or lesser extent by each of the following: the direction of greatest sub-laminate stiffness, global buckling/bending and the direction of the fibres in the bottom ply of the sub-laminate. If the base laminate were to remain completely flat, it is likely that the direction of propagation would be mostly influenced by the maximum stiffness direction of the sub-laminate and bottom ply fibre direction. However, most laminates displayed global buckling/bending, which encouraged propagation to occur transverse to the load (figure 8), and some laminates were influenced by the direction of the fibres in the bottom ply of the sub-laminate, which encouraged growth parallel to the bottom ply fibre direction. This is because growth perpendicular to the fibre direction was inhibited by intra-ply cracking. Figure 5 shows how the three mechanisms combine to affect the direction of propagation in the [45_{2}] laminate. Initiation of propagation occurs at a point on the perimeter of the PTFE insert away from the direction of the greatest sub-laminate stiffness at an angle of approximately 60° (figure 5*c*) owing to the influence of global buckling/bending and then subsequently propagates in the 45° direction as intra-ply cracking prevents growth in any other direction.

Intra-ply cracking was noted in sub-laminates having two plies of identical orientation on their upper surface (figure 5). In these cases, the intra-ply cracks drove the direction of growth as they prevented delamination propagation occurring perpendicular to the fibre direction. In the [90/0/90] and [0_{2}/90] sub-laminates, delamination growth occurred at the artificially delaminated interface but it was not confined to this interface. Instead, following initial propagation at the artificially delaminated interface, intra-ply cracking allowed growth to transfer to the second ply interface in the sub-laminate. Although this is an interesting phenomenon, it occurred following initial propagation and so is outside the capabilities of the strip and FEA models and hence is not discussed further here. Similar mechanisms have been recorded by Greenhalgh and Singh [17].

### (c) Finite element results

The parameters used for cohesive elements within the FEAs are summarized in table 3 and table 2 gives buckling, initiation and propagation results. The values for stiffness and mode I strength are from Lachaud *et al*. [18] and the mode II and III strengths were supplied directly by the Hexcel Corporation. The value for *G*_{IIC} was chosen to represent a lower bound value based on initial fracture of a 45° plane under pure shear . This value was subsequently confirmed to be similar to the value obtained by experimental characterization [19]. The value of *η* was chosen to represent a significant mode II contribution (see equation (2.20)) and is therefore lower than some other values used in the literature [20]. As an example of local buckling output from Abaqus, figure 9 shows the critical buckling mode of the [0/90_{2}] laminate obtained from linear eigenvalue analysis. The delamination area buckles with two half-wavelengths, which agrees with the experimental buckling mode shown in figure 7*a*. However, owing to contact with the base laminate, the experimental sub-laminate buckle forms with unequal amplitude for each of the half-waves in the load direction; this contact problem is not captured by the linear FEA model as superposition of the base laminate and sub-laminate was not prevented. However, contact is modelled during the nonlinear propagation analysis presented below.

For the [45_{2}] laminate, propagation initiates at about 60° to the load direction, as shown in figure 10*a*. This is coincident with the experimental test (see figure 5*c*), although intra-ply cracking subsequently forced experimental propagation to occur in the 45° direction (figure 5*d*). Figure 10*b* shows the delamination propagation of the [0_{2}/90] laminate simulated by the Abaqus model. As noted in figure 10, Abaqus predicts that delamination initiation occurs when the nodes start to separate at 2838 μstrain. Propagation is predicted to occur transversely along a horizontal line through the centre of the delamination starting from a point on the perimeter of the delamination. Increasing applied strain in the Abaqus model leads to propagation where a delamination area is formed as full elements fail (each element has four nodes, which are all required to fail for the element to fail), as depicted in figure 10. The delamination then propagates in the horizontal direction, which is shown by the failure of cohesive elements.

With the exception of the [45_{2}], [45/−45] and [0/45] laminates, where initiation strains were significantly unconservative, both the delamination initiation strain and subsequent directions of delamination propagation predicted by Abaqus are in good agreement with the experimental tests. As shown in table 2, propagation strains are between 5 and 65 per cent higher than initiation strains. For the cross-ply laminates [02], [90/0/90], [0/90_{2}], [90_{2}/0] and [0_{2}/90], propagation strains are between 16 and 37 per cent higher, whereas for the angle-ply laminates [45_{2}] and [45/−45], they are approximately 5 per cent higher. The largest difference of 65 per cent was noted for the most asymmetric sub-laminate, [0/45]. This divergence indicates that the stacking sequence of the sub-laminate affects the mode mixity. The parameter *η* in equation (2.20) is defined to model the mode mixity and should be derived from curve fitting to experimental data obtained from tests defined to capture pure modes I and II, and mixed-mode behaviour [20]. This implies that the value of *η*=1.75 used here may not be suitable in all cases.

### (d) Strip model results

Strip model local buckling and propagation threshold strain results are also given in table 2. Figure 11 shows a local buckling mode shape produced by Viconopt [10] for the [0/90_{2}] sub-laminate which corresponds to the equivalent experimental (figure 7) and FEA (figure 9) local buckling mode shapes. Viconopt does not take into account contact between the sub-laminate and base and hence does not produce the unequal half-wave amplitudes seen in the experiments, such as the [0/90_{2}] experiment (note the number of contours on each half-wave of the buckle in figure 7*a*), where sub-laminate and base contact occurs. However, a comparison of strip model and experimental results in table 2 suggests that contact between the sub-laminate and base during local buckling has limited influence on the accuracy of the buckling strain predictions. In fact, where the critical buckling strain *ε*^{C}>*ε*_{th}/3 the use of a lower *ε*^{C}, violating contact is always conservative compared with a higher, non-contacting mode. For example, for a [0/90_{2}] sub-laminate the first Viconopt buckling mode having one full wave (i.e. without violating the contact condition) occurred at 2320 μstrain.

Use of this *ε*^{C} gives *ε*_{th}=3780 μstrain, a 2–3% increase compared with the value given in table 2.

In general, propagation strains produced by the model are in agreement with both experimental and FEA results. The only exceptions to this being the [0_{2}/90] and [0/45] sub-laminates which display, respectively, local–global interaction and full elastic coupling, which are not captured by the strip model, and the [45_{2}] and [45/−45] laminates for which the strip model gives considerably more accurate predictions of experimental propagation initiation than does the FEA. However, for the [0_{2}/90] case, as noted in §4, local strains detected by the DIC in the areas into which propagation subsequently occurred were equal to the strains predicted by the strip model.

## 4. Discussion

In the case of the [45_{2}], [0/90_{2}], [0_{2}/90] and [0/45] sub-laminates experimental buckling occurred without adhesion and hence all buckling strain results are well correlated. The [0_{2}], [45/−45], [90/0/90] and [90_{2}/0] sub-laminates suffered from residual adhesion between the base laminate and sub-laminate; the analytical models do not account for this and hence in these cases local buckling results are not well correlated with those from the experiments. Such adhesion did not seem to affect accuracy of either the FEA or the strip model in predicting threshold strains. This is because, provided the PTFE region is eventually released, similar amounts of strain energy are developed. Other phenomena such as (i) contact between sub-laminate and base during local buckling, (ii) intra-ply cracking, and (iii) propagation at multiple interfaces did not seem to affect the accuracy of the model, though in the last case this is not unexpected as growth initiates at a single interface and does not spread to other interfaces until this has occurred. Note in particular that surface ply cracking the [0_{2}], [45_{2}] and [90_{2}/0] laminates did not appear to affect propagation strain predictions or growth direction as there was minimal stiffness and energy stored transverse to the fibre direction in these sub-laminates.

As in §3*b*, some coupons demonstrated global bending. For example, figure 6 displays global buckling as a large change in gradient of the strain gauge curves between 100 and 140 kN, indicating a change in stiffness. For most cases in table 2, interaction between local and global buckling modes has little effect on threshold strain. However, this is not true of the [0_{2}/90] coupon, which has a considerable amount of axial stiffness in its delaminated sub-laminate. Following buckling of the sub-laminate, the bending stiffness of the [0_{2}/90] laminate was considerably reduced and hence bending of the base occurred at a level of applied strain that was significantly lower than for other laminates. This is demonstrated by comparison of figure 5*b* showing limited global deformation (wide outer contours) of the [45_{2}] laminate and figure 8*b* showing significant global deformation (tight outer contours) in the [0_{2}/90] laminate at an equal level of load. The large out-of-plane deformations displayed by the [0_{2}/90] caused a considerable local increase in compressive strain around the perimeter of the local buckle (figure 12) that was not detected in the (far-field) strain measurements collected by the strain gauges. Figure 12 shows that the local strain in the region into which the delamination grows at initiation is approximately equal to the strip model prediction. In the case of the [0_{2}/90] laminate a large area has local strain greater than the threshold strain and hence propagation continues at a significant rate with increasing applied displacement. However, if the area into which propagation could occur is small (i.e. a very localized area with strain above threshold for the laminate) then growth will not occur as swiftly. Hence interaction between local and global buckling can have an effect on the rate at which propagation occurs and in cases where the areas of local strain above the threshold are very localized they are unlikely to be critical provided the far-field strain remains below the threshold strain. An example of this is the [0/45] laminate, although in this case the propagation is also likely to have been influenced by the local post-buckling mode shape of the fully coupled sub-laminate which displays bend–twist coupling and is both asymmetric and unbalanced.

Note that for the purpose of optimum laminate design the results of table 2 indicate that delaminations creating sub-laminates containing ±45° angle plies are damage tolerant—in other words, they display the greatest strength–whereas delaminated sub-laminates containing 0° ply are not damage tolerant. These observations are easily apparent from application of the model in the context of optimum design.

The initiation of damage predicted by FEA is generally closer to the experimental propagation values than the FEA propagation strains. Use of the latter would give a prediction of propagation that is up to 85 per cent higher than the experimental value. For all but the [0_{2}/90] and [0/45] cases considered above, the strip model actually gives better agreement with experimental values than FEA initiation. This may be due to the inadequacy of the mixed-mode formulation of equation (2.20) to capture all combinations of mode mixity.

## 5. Conclusions

New developments of a simple and efficient strip model for predicting the threshold propagation strain of compressively loaded, delaminated composites have been presented. The accuracy and efficiency of the model make it well suited to application in optimum ply-stacking algorithms designed to maximize laminate strength. Indeed, a comparison of results with those from FEA and experimental compression tests indicates that, in most cases, the strip model can conservatively predict the strain at which delamination propagation occurs to within 10 per cent of experimental values and 11 per cent of FEA initiation values. The quoted accuracy of the model is true for sub-laminates displaying a range of anisotropy, including 90° ply dominated sub-laminates, subject to the following conditions: (i) the thin-film assumption made in the modelling methodology must hold and (ii) full elastic coupling effects must not play a significant role in the post-buckling of the sub-laminate. Nonlinear FEA prediction of propagation using cohesive elements was shown to be significantly inaccurate in some cases, although initiation values gave reasonable agreement.

Future work should consider the implications of full elastic coupling in the sub-laminate and the interaction between local and global buckling, particularly if global buckling occurs before initial propagation.

## Acknowledgements

The authors are grateful for discussions with Prof. Giles W. Hunt and David Lovell during preparation of the manuscript. The second and third authors are currently supported by the EPSRC (EP/H025898/1), Airbus Operations and GKN Aerospace.

## Footnotes

One contribution of 15 to a Theme Issue ‘Geometry and mechanics of layered structures and materials’.

- This journal is © 2012 The Royal Society