## Abstract

A damage-based finite-element model is used to predict the fracture behaviour of centre-notched quasi-isotropic carbon-fibre-reinforced-polymer laminates under multi-axial loading. Damage within each ply is associated with fibre tension, fibre compression, matrix tension and matrix compression. Inter-ply delamination is modelled by cohesive interfaces using a traction-separation law. Failure envelopes for a notch and a circular hole are predicted for in-plane multi-axial loading and are in good agreement with the observed failure envelopes from a parallel experimental study. The ply-by-ply (and inter-ply) damage evolution and the critical mechanisms of ultimate failure also agree with the observed damage evolution. It is demonstrated that accurate predictions of notched compressive strength are obtained upon employing the band broadening stress for microbuckling, highlighting the importance of this damage mode in compression.

This article is part of the themed issue ‘Multiscale modelling of the structural integrity of composite materials’.

## 1. Introduction

Over the last four decades, significant efforts have been made to develop predictive tools for composite laminate failure, and a significant literature now exists on the tensile and compressive responses. For example, Awerbuch & Madhukar [1] reviewed early fracture models (such as the point stress criterion) for the tensile strength of composite laminates containing stress-raisers. Subsequent models included the role of subcritical damage (such as splits and delamination) upon tensile strength (e.g. [2–4]). By quantifying the blunting effect of notch tip damage, the tensile strength of cross-ply carbon-fibre-reinforced-polymer (CFRP) specimens was predicted for a range of notch sizes and lay-ups. Finite-element (FE) simulations followed. For example, Wisnom & Chang [4] used interface/cohesive elements to model inter-ply delamination and intra-ply splits in cross-ply tensile specimens. Hallett *et al*. [5] extended this model using a Weibull-based fibre failure criterion to predict the strength of notched quasi-isotropic CFRP specimens. Parallel efforts to model tensile fracture include the micro-scale approach [6], meso-scale approach [7] and multi-scale strategy [8,9]. More recently, user-defined constitutive models have been used to predict progressive tensile failure (e.g. [10]).

Damage models have also been generated for compressive failure [11–14]. For example, Soutis and co-workers [11] have developed a cohesive zone model wherein notch tip damage (0° microbuckling, delaminations, etc.) is represented by an equivalent crack loaded on its faces; their model was used successfully to predict the effects of hole size and lay-up on the compressive strength of a CFRP. More detailed approaches followed, based on fibre kinking [12–15]. By contrast, a limited literature exists on the measurement and prediction of shear failure. Hollmann [16] developed a damage zone model and obtained accurate predictions of shear strength, but did not attempt to predict the damage modes. There remains a clear need to develop an experimentally validated tool to predict damage development and failure under multi-axial loading, thereby motivating the present paper.

In this study, the Hashin–Rotem damage model is used to predict the progressive damage and failure of quasi-isotropic CFRP laminates, containing either a notch or a circular hole and subjected to multi-axial loading. The intent is to compare predictions against recent experimental observations [17] for the geometries shown in figure 1. Four types of damage are modelled: fibre tension, fibre compression, matrix tension and matrix compression. Inter-ply delamination is represented by cohesive interfaces and an assumed traction-separation law. Predictions are compared with the recent data [17], with an emphasis upon (i) failure envelopes under combined tension and shear, (ii) ply-by-ply and inter-ply subcritical damage evolution, and (iii) the mechanisms that dictate ultimate failure. Three distinct failure mechanisms were observed in [17] and are summarized in figure 2. The observed mechanism depends upon the loading state (combined tension, compression and shear), and in each case the stress concentration is reduced by limited splitting at the notch root:

(i) Mechanism A: the application of direct tension, or combined tension and moderate shear, leads to tensile rupture of the main loading-bearing 0° plies, with splitting of the 0° plies.

(ii) Mechanism B: under pure shear, or shear and moderate tension, microbuckling of the main loading-bearing −45° plies occurs, with splitting of the −45° plies.

(iii) Mechanism C: under direct compression, microbuckling of the main loading-bearing 0° plies occurs, with splitting of the 0° plies.

The challenge is to predict each of these failure mechanisms and their domain of dominance.

## 2. Finite-element model

Three-dimensional FE simulations of the progressive damage and failure of centre-notched quasi-isotropic CFRP laminate specimens under multi-axial loading were performed using the explicit version of the commercial FE package Abaqus FEA (v. 6.12). The specimen geometry was identical to that in [17], and the gauge section is sketched in figure 1. Quasi-isotropic [45/0/−45/90]_{2S} laminate specimens, of thickness 2 mm and of rectangular gauge section 23×25 mm, were made from HexPly^{r} IM7/8552 CFRP. Two types of centre notch were assumed: (i) a circular hole of diameter 6 mm or (ii) a notch of length *a*=6.25 mm, height *b*=0.70 mm and root radius *b*/2. Displacements were applied to the periphery of the mesh to generate combined tension and shear, and compression and shear. Preliminary calculations reveal that each specimen responded with proportional remote stressing up to the point of catastrophic failure. Define the average shear traction by *τ* and the tensile normal traction by *σ* as done in [17]. Then, the loading direction *φ* is defined via tan *φ*=*τ*/*σ*.

Each ply was modelled using a single layer of 8-noded reduced-integration continuum shell elements (SC8R), and delamination at the interface between adjacent plies was modelled using a single layer of 1 μm thick cohesive elements (COH3D8). The ‘general contact’ option in Abaqus was employed to simulate contact between all possible surfaces of the plies. The mesh density was greatest at the stress raiser in order to capture the elastic stress concentration factor *k*_{T}; to achieve this, the smallest element was of side length 50 μm.

### (a) Constitutive model for each ply

The anisotropic damage model within Abaqus Explicit is based on the work of Hashin & Rotem [18], Hashin [19], Matzenmiller *et al*. [20] and Camanho & Davila [21]. The initially undamaged material exhibits a linearly elastic response. Hashin’s initiation criteria are used to predict the onset of damage, and a fracture energy-based constitutive law with linear material softening is used to predict damage evolution. Four different modes of failure are accounted for: (a) fibre failure in tension, (b) fibre microbuckling in compression, (c) matrix cracking under transverse tension and shear, and (d) matrix crushing under transverse compression and shear.

Consider a fibre-reinforced unidirectional ply, with the 1-axis along the fibre direction, the 2-axis in the in-plane transverse direction and the 3-axis normal to the plane of the ply. The ply is transversely isotropic with respect to the fibre direction. The Young’s moduli *E*_{1} and *E*_{2} in the 1- and 2-directions, along with the in-plane shear modulus *G* and Poisson ratio *v*_{21} are the relevant elastic constants of the unidirectional lamina.

The laminate is taken to be linear elastic, as specified above, up to the initiation of damage. A nonlinear stress versus strain response accompanies damage progression due to a progressive drop in the three moduli (*E*_{1}, *E*_{2}, *G*) with increasing strain. Four scalar damage variables are introduced, corresponding to the four damage modes (a)–(d) as listed above. In the undamaged state, each damage variable is set to zero. As strain-controlled damage accumulates, one or more damage variables increase to a maximum value of unity. The moduli drop with increasing value of damage variables, such that one or more moduli vanish when one of the damage variables attains unity. The damage evolution law follows that laid down by Matzenmiller *et al*. [20].

The damage evolution law comprises two steps. First, damage initiates when a critical strain state is attained, as proposed by Hashin [19]. This condition is usually re-expressed in terms of the associated critical stress state akin to the yield surface in plasticity theory: for a stress state within the damage locus, no additional damage occurs, and the stress versus strain response is linear and reversible. Second, the damage variable(s) increase nonlinearly with increasing strain, and this leads to a drop in moduli and thereby to a drop in stress.

Denote the tensile and compressive strengths for damage initiation in the fibre direction (1-direction) by *X*^{T} and *X*^{C}, respectively. The tensile and compressive strengths in the transverse direction are denoted by *Y* . After damage has developed, these strengths drop as follows. Write the damage variable for tensile failure in the fibre direction as *d*^{t}_{f}. Then, the current tensile strength in the fibre direction equals Likewise, the damage variable for compressive loading in the fibre direction is *d*^{c}_{f}, while those for transverse tension and compression are *d*^{t}_{m} and *d*^{c}_{m}, respectively. The damage remains constant when the stress state lies within the following critical surfaces:
2.1*a*
2.1*b*
2.1*c*
2.1*d*where 〈〉 represents the Macaulay bracket, such that 〈*x*〉=0 for *x*<0 and 〈*x*〉=*x* for *x*≥0. The shear damage *d* is defined as
2.2In any given state of damage, the secant relationship between stress and strain reads
2.3

where
2.4*a*and
2.4*b*Note that the composite can sustain tensile fibre stress when fully damaged in, say, the compressive mode, and so on. The validity of this assumption remains to be verified experimentally. The evolution law is summarized in appendix A for each of the four independent damage variables , d^{c}_{f}, *d*^{t}_{m} and *d*^{c}_{m}. The damage growth law adopted here is based upon the assumption that the stress decreases linearly with increasing strain once damage initiates. Table 1 lists the material parameters inputted into the model for each composite ply. They are all experimentally measured values collated from the literature on the HexPly^{r} IM7/8552 material system.

### (b) Cohesive law for the interface between plies

The three-dimensional 8-noded cohesive element (COH3D8) was used to attach the unidirectional plies in the laminates to each other and its traction-separation behaviour allowed the FE model to capture delamination between plies. Define the three components of traction in the cohesive layer by the normal traction *t*_{n} (along the local 3-direction) and the two shear tractions, *t*_{s} (along the local 1-direction) and *t*_{t} (along the local 2-direction). The corresponding normal and shear separations are denoted by *d*_{n}, *d*_{s} and *d*_{t}. The traction-separation relation of the cohesive element is of the form
2.5*a*
2.5*b*
2.5*c*where is the normal stiffness in terms of the peak normal traction *t*^{0}_{n} and the corresponding normal displacement *δ*^{0}_{n}. The corresponding shear stiffnesses are analogously defined as and The scalar damage variable *D* lies in the range 0≤*D*≤1 and represents the overall damage in the material. *D* equals zero for the initial, undamaged state and *D* equals unity for the fully damaged state. *D* increases linearly with increasing effective cohesive displacement *δ*_{m} as defined by
2.6where 〈*x*〉=*x* for *x*≥0 and 0 otherwise. It remains to define the state at *D*=1.

The ‘energy release rate *G*_{n} in the normal direction’ is defined by the area under the normal traction–displacement curve to the current state. The shear components *G*_{s} and *G*_{t} are defined likewise. At final fracture (*D*=1), the following relationship is satisfied:
2.7in terms of the measured toughness in the normal mode and in the two shear modes *G*^{C}_{s} and *G*^{C}_{t}. This is consistent with the data of Li *et al*. [29].

Table 2 lists the input parameters for the cohesive elements. With the exception of the stiffnesses, all input parameters are taken from measurements in the literature for the HexPly^{r} IM7/8552 system. The choice for the stiffnesses of the cohesive layer is arbitrary but is sufficiently large for numerical stability. This was achieved by taking *δ*^{0}_{n} to be 0.1 μm, and *δ*^{0}_{s} and *δ*^{0}_{t} to be 1 μm.

## 3. Results

The FE simulations were used to predict the failure envelope, and the ply-by-ply and inter-ply subcritical damage evolution for the notch and circular hole geometries of figure 1. These predictions were compared with the detailed experimental data of Tan *et al*. [17] in order to investigate the reliability and accuracy of the numerical approach. Having validated the FE model, additional predictions were performed in the compression-shear regime, for which experimental data are more limited.

### (a) Failure envelope

The predicted failure of quasi-isotropic CFRP laminate specimens containing a central notch or hole was accompanied by a dynamic drop in load. The peak strength was plotted in tensile/compressive-shear stress space to construct the failure envelopes; see figure 3*a*,*b* for the notch and hole, respectively.

We emphasize that there is choice in the most appropriate value of longitudinal compressive strength *X*^{C} of the plies to be used in the FE model. Recall that geometric imperfections cause the (nominal) peak compressive stress associated with microbuckling to drop from a peak strength to the steady-state band broadening stress (e.g. [25]). It is hypothesized here that the band broadening stress is the most suitable value for *X*^{C}. To test this, a compression test on a specimen with a notch was repeated by first sharpening the notch root by a razor blade. It was found that the compressive strength of the specimen with the sharpened notch equals that of a specimen with an unsharpened notch. Further, the length of microbuckle prior to ultimate failure, as observed through X-ray computed tomography (CT), was the same for both specimens. This suggests that the as-manufactured central notch induces a sufficiently large imperfection in terms of fibre waviness that the peak stress approaches the band broadening stress. Consequently, the band broadening stress value for *X*^{C} was used in our model. According to Moran *et al*. [25], the post critical steady-state band broadening stress of a composite with an effective shear strength of 100 MPa can range between 200 and 800 MPa. An average value within that range was employed *X*^{C}=500 MPa in this study.

The significance of the assumed value of compressive strength *X*^{C} is explored in figure 3*a* for the case of a notch. Simulations were performed by taking *X*^{C}=1730 MPa (unnotched compressive strength) and *X*^{C}=500 MPa (band broadening stress). It is clear from figure 3*a* that the more accurate predictions are obtained by assuming the compressive strength is set by the band broadening strength *X*^{C}=500 MPa. Consequently, for the remainder of our study, we assume that *X*^{C}=500 MPa for both the notch and circular hole configurations.

In broad terms, predictions for the failure envelope of notch and circular specimens are in good agreement with the measurements (figure 3). However, there are a few minor discrepancies:

(i) The FE model predicts the angle

*φ*at which Mechanism A transitions to Mechanism B for the circular hole configuration equals 45°. By contrast, experiments reveal a transition angle of*φ*=30° (figure 3*b*).(ii) The model predicts a slightly larger degree of interaction between Mechanism A and Mechanism B than experimentally observed for the notch specimen, as evidenced by the curved profile of the predicted failure envelope in figure 3

*a*.

#### (i) Predictions for compression-shear

Jelf & Fleck [31] observed that unidirectional composites fail by plastic microbuckling under combined compression and shear; the axial compressive strength decreases linearly to zero as the longitudinal shear stress is increased to the composite shear strength. Edgren *et al*. [32] observed a similar response for quasi-isotropic non-crimp fabric laminates under combined compression and shear. The FE model of this study predicts that the open-hole and notch specimens fail by mechanism C under combined compression and shear. And for both geometries, the compressive strength decreases in an almost linear fashion with increasing shear stress, as witnessed by the near-linear slope of the failure envelopes in figure 3.

### (b) Damage evolution: prediction versus observation for notch specimens

It is instructive to compare the predicted ply-by-ply evolution of damage with the experimental observations of Tan *et al*. [17]. The qualitative nature of the observed damage is the same for the notch and hole (with slightly greater damage development for the hole associated with its larger stress concentration factor), and so it suffices to give a detailed comparison here of predicted versus observed damage development for the notch. We consider damage development by the overall effective measure *d* (as defined in (2.2)) for mechanisms A, B and C in turn, and then report the role of delamination for all three mechanisms. A brief summary of damage development in the hole geometry is then presented to contrast it with the notch case.

#### (i) Mechanism A

The experiments of Tan *et al*. [17] reveal the presence of extensive intra-ply splitting in all plies of the laminate for the case of remote tensile loading (*φ*=0°). Excellent correlation is achieved between the predicted and measured subcritical damage modes in each ply of the quasi-isotropic laminate, as shown in figure 4*a*. The 0° splits are significant as they reduce the stress concentration of the notch and hole, and delay the onset of tensile rupture of the 0° fibres. As seen in figure 4*a*, the extent of these 0° splits is somewhat underestimated by the FE model, but despite this, the failure strength of the specimen at *φ*=0° is accurately predicted, recall figure 3*a*.

The FE model reveals that intra-ply splitting in all plies initiates at approximately 18% of the failure load. However, the splits in each ply orientation grow at different rates: the splits in the 90° plies grow the fastest, followed by the ±45° splits and then the 0° splits. This explains why the 90° splits are first detected experimentally, followed by the ±45° splits and then the 0° splits [17] (the X-ray CT technique has a finite resolution). The limited amount of fibre tensile fracture that occurs prior to ultimate failure in the load-bearing 0° plies is also captured by the model (figure 4*b*).

When *φ* equals 15°, the loading is predominantly tensile with a small level of shear stress. As discussed by Tan *et al*. [17], a limited degree of superimposed shear stress causes the splits in the load-bearing 0° plies to become asymmetrical. This feature is captured well by the FE model (figure 5) which gives a comparison of the predicted and measured subcritical damage modes at 90% failure load. The critical failure mechanism for Mechanism A is fibre tensile failure of the load-bearing 0° plies.

Figure 6 shows a comparison of the predicted 0° split length versus load with the measured response at *φ*=15°. A good correlation is evident. The splits initiate very early in the loading process, and once initiated, they grow approximately linearly with increasing remote load.

#### (ii) Mechanism B

Recall from Tan *et al*. [17] that the dominant damage mechanism associated with Mechanism B is microbuckling from the notch tips in the load-bearing −45° plies. We now demonstrate that the FE model accurately captures this feature and also predicts its evolution to acceptable accuracy. For example, the extent of microbuckling in the −45° plies of the notch specimen under pure shear appear is adequately predicted by the model (figure 7). However, the direction in which the microbuckles grow deviates from the observed direction. In the experiments, the microbuckles advance in a straight trajectory directly ahead of the notch, whereas the simulations predict that the microbuckles grow at an inclination to the notch. The reason for this discrepancy is not immediately clear. One possible explanation is that the model does not fully capture the effect of adjacent splits upon the direction of microbuckle propagation. As suggested by Tan *et al*. [17], splitting in the adjacent plies plays a significant role in guiding the direction of microbuckle propagation.

Excellent correlation of predicted and measured microbuckle length is also noted for *φ*=30° (figure 8). In summary, the FE model correctly predicts that the critical failure mechanism for Mechanism B is plastic microbuckling of the load-bearing −45° plies. Failure ensues when a microbuckle grows across the ligament of the specimen.

#### (iii) Mechanism C

Recall from Tan *et al*. [17] that the dominant damage mechanism for Mechanism C is plastic microbuckling of the load-bearing 0° fibres. The FE model accurately predicts this mode of damage and its evolution with increasing remote load (figure 9). The predicted microbuckle lengths correlate well with the observed microbuckle lengths at each level of loading. The model also correctly predicts that ultimate failure occurs when the microbuckles in the 0° plies propagate across the ligament of the specimen.

#### (iv) Delamination

Interfacial delamination for each of the three mechanisms is predicted by the model via the cohesive elements. The observed state of delamination at about 90% of failure load has been measured by X-ray CT by Tan *et al*. [17] and is compared with predictions in figure 10. In broad terms, the FE simulations adequately capture the delamination state. For Mechanism A, the observed and predicted delamination zones are localized at the notch tips and bridge adjacent intra-ply splits. For Mechanism B, the delamination is more extensive and occurs mostly around the microbuckled region of the −45° plies. Finally, for Mechanism C, delamination accompanies the limited amounts of 0° fibre microbuckling prior to ultimate failure.

### (c) Damage evolution: prediction versus observation for hole specimens

The predicted and observed damage mechanisms are compared in figure 11 for the critical ply of each mechanism, at 90% failure load for the circular hole: good correlation is evident. For Mechanism A, the observed (and predicted) splits was short at the edge of the hole in the 0° ply orientation. For Mechanism B, the splits and microbuckles that develop at the edge of the hole in the −45° plies are adequately predicted. Likewise, the microbuckle in the 0° plies for Mechanism C is predicted by the model.

Now consider a representative case of the ply-by-ply damage state in a circular hole specimen. A comparison of the predicted and measured ply-by-ply damage in the circular hole specimen at *φ*=30°, and loaded to 90% failure load (Mechanism B), is given in figure 12. It is clear that the FE model gives an accurate prediction for the damage state within each of the plies.

## 4. Concluding discussion

In the literature, the prediction of progressive damage and failure in notched quasi-isotropic laminates is largely focused upon the pure tension case [5,10,33]. A brief comparison between the accuracy of our model and alternative models in this respect is now discussed. Hallett *et al*. [5] and Li *et al*. [33] developed an FE approach which explicitly modelled delamination and splits using interface elements, and fibre failure using a progressive statistical failure theory. This approach assumed the form and location of damage *a priori*. Owing to the detailed modelling of the individual damage modes, the prediction of tensile damage was accurate. The extent of the intra-ply splits, including the 0° splits, was captured with high accuracy. However, the accurate FE strategy of Hallett *et al*. [5] and Li *et al*. [33] is traded against the simplicity of our model for the wide range of subcritical damage mechanisms exhibited in Mechanisms A, B and C. In particular, the failure strength of the tensile specimen is adequately predicted by our model (recall figure 3). This suggests that the current approach is sufficient to yield accurate strength predictions.

Chen *et al*. [10] developed a combined elastoplastic damage model which accounts for both plasticity and damage to represent the mechanical response of composite layers through a user-defined material subroutine. Delamination was accounted for by using cohesive layers. Chen *et al*. [10] found excellent correlation with experiments for the progressive failure of [0/45/90/−45]_{2S} and [0_{2}/45_{2}/90_{2}/−45_{2}]_{S} AS4/PEEK laminates with central circular holes subjected to in-plane tension. They compared their elastoplastic model with the Abaqus in-built progressive damage model (as employed herein) and found that the analyses using the Abaqus model overestimate the tensile failure loads. However, Chen *et al*. [10] took the unnotched compressive strength of the ply to be the longitudinal compressive strength *X*^{C} parameter. Recall that an overestimation of the tensile strength was also obtained in the present study when the unnotched compressive strength was adopted for *X*^{C} (figure 3). However, upon switching the unnotched compressive strength to the band broadening stress for *X*^{C}, we obtained excellent agreement between the predicted and measured open-hole (and notched) tensile strength of the laminate.

It is recognized that small discrepancies exist between the simulations of the present study and the measurements of Tan *et al*. [17], and this is due to several factors:

(i) The assumed value for the band broadening stress was estimated and not measured.

(ii) The current constitutive law assumes linear softening to the plies, once damage has initiated. This is idealized and may not be representative of the actual strain softening response.

(iii) Laminate ply-by-ply failure analysis is performed on the basis of the two-dimensional stress field in the laminate. In reality, more complex damage initiation criteria could be employed such as the LaRC02 criteria developed by Davila

*et al*. [34], an enhanced version of Hashin’s failure criteria.(iv) The FE model lacks some of the sophistication of other numerical approaches such as in Hallett

*et al*. [5], Li*et al*. [33], Daghia & Ladevèze [8] and Chen*et al*. [10]. Nevertheless, the model does capture the main damage mechanisms over a wide range of loading states.

This study is in broad support of a simple FE approach that uses experimentally ascertained material parameters as inputs. The model has some success in predicting the complex multi-axial fracture behaviour of notched quasi-isotropic CFRP laminates. Its ability to predict the role of lay-up in multi-axial strength will be explored in a subsequent paper.

## Authors' contributions

J.L.Y.T. carried out the laboratory work and participated in data analysis; V.S.D. participated in the design of the study and drafted the manuscript; N.A.F. conceived of the study, designed the study, coordinated the study and helped draft the manuscript. All authors gave final approval for publication.

## Competing interests

The authors declare that they have no competing interests.

## Funding

We received no funding for this study.

## Appendix A. The damage evolution law

First, consider the tensile fibre damage mode. An effective strain is defined by and is used to update the damage state variable via the relation
A 1where X^{T}/*E*_{1} is the value of when the tensile fibre first initiates and 2*J*^{t}_{f}/*L*_{e}X^{T} is the strain for a complete tensile fibre damage. Here, *J*^{t}_{f} is the tensile fibre fracture energy, and *L*_{e} is a representative length scale. When the axial strain is decreased, we hold *d*^{t}_{f} fixed. Similarly, the compressive fibre damage variable is specified as
A 2where *J*^{c}_{f} is the compressive fibre fracture energy and the effective strain is The matrix damage variables are given by
A 3and
A 4where *J*_{m} is the matrix fracture energy and the effective strains are defined as and In this study, we choose *L*_{e} to be the size of an FE; numerical experimentation confirmed that this choice gives a response that is almost independent of mesh size.

After the failure criterion for a particular mode is met, the damage variable evolves according to the following equation: A 5The damage variable does not exceed a value of unity. The model is designed to remove or delete an element once the damage variables for all failure modes at all material points reach unity.

## Footnotes

One contribution of 22 to a Theo Murphy meeting issue ‘Multiscale modelling of the structural integrity of composite materials’.

- Accepted February 24, 2016.

- © 2016 The Author(s)

Published by the Royal Society. All rights reserved.