## Abstract

The paper presents an analytical approach to predicting the effect of intra- and interlaminar cracking on residual stiffness properties of the laminate, which can be used in the post-initial failure analysis, taking full account of damage mode interaction. The approach is based on a two-dimensional shear lag stress analysis and the equivalent constraint model of the laminate with multiple damaged plies. The application of the approach to predicting degraded stiffness properties of multidirectional laminates under multi-axial loading is demonstrated on cross-ply glass/epoxy and carbon/epoxy laminates with transverse and longitudinal matrix cracks and crack-induced transverse and longitudinal delaminations.

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

## 1. Introduction

Failure process of fibre-reinforced composite laminates subjected to multi-axial loading involves sequential accumulation of intra- and interlaminar damage in the form of matrix cracking and delamination. Intralaminar matrix cracks parallel to the fibres in the off-axis plies is the first damage mode observed. Depending on the laminate stacking sequence, these cracks are either arrested at the interface or cause interfacial failure leading to delamination and/or cracking in the adjacent layers due to high interlaminar stresses at the interface.

Development of intra- and interlaminar damage in composite laminates has been the subject of numerous studies in the literature (e.g. [1,2]). More recently, Rebiere & Gamby [3,4] proposed an energy criterion based on the computation of the partial strain energy release rates associated with transverse cracking, longitudinal cracking and crack-induced delamination and used it to predict initiation of these damage modes in symmetric cross-ply laminates subjected to uniaxial loading. Lim & Li [5] evaluated the energy release rates associated with matrix cracking and crack-induced delamination and used them to critically evaluate damage mode transition from transverse cracking to delamination. Blazques *et al.* [6] employed boundary element method to carry out a numerical study of the stress state in the neighbourhood of matrix crack-induced delamination in a cross-ply laminate in order to clarify the mechanisms of damage interaction between transverse cracking and delamination. Maimi *et al.* [7,8] carried out a comprehensive study of matrix cracking and crack-induced delaminatons, proposing a model to simulate stress–strain state of the damaged ply and using it to analyse evolution of matrix cracking and crack-induced delamination. Garcia *et al.* [9] modelled transverse cracking onset and growth in cross-ply laminates using a coupled stress and energy criterion. Intra- and interlaminar cracking in composite laminates under impact loading and four point bending was investigated by Shi *et al.* [10] and Shi & Soutis [11], respectively.

Multidirectional laminates subjected to multi-axial loading may still be capable of carrying load after matrix cracking has occurred. In the laminate, in-plane shear and normal stresses can be transferred, to some extent, back into the damaged lamina via the neighbouring laminae. Owing to this stress transfer, damaged lamina within the laminate retains a certain amount of load-carrying capacity. *In situ* stiffness of a damaged lamina constrained within the laminate depends on the damage configuration and stiffnesses and thicknesses of neighbouring laminae. Prediction of the post-initial failure behaviour of a laminate requires accurate information regarding the properties of the damaged lamina.

This paper describes a method of predicting the effect of intra- and interlaminar damage on the stiffness properties of the laminate which can be used in the post-initial failure analysis, taking full account of damage mode interaction. The approach is based on the equivalent constraint model (ECM) of the damaged laminate [12–21]. Closed-form expressions are given for the *in situ* damage effective functions which characterize degraded stiffness properties of each damaged ply; for a given damaged ply they explicitly depend on the damage parameters (matrix crack density and relative delamination area) associated with that ply and implicitly on the damage parameters associated with other damaged plies.

The application of the approach to predict the degraded stiffness properties of multidirectional laminate with multilayer inter- and interlaminar damage is shown for cross-ply glass/epoxy and carbon/epoxy laminates damaged by transverse and longitudinal matrix cracks and crack-induced transverse and longitudinal delaminations.

## 2. Equivalent constraint model

Figure 1 shows a schematic of the cross-ply [0_{m}/90_{n}]_{s} laminate damaged by transverse and longitudinal delaminations growing from the tips of transverse cracks in the 90° plies and longitudinal cracks in the 0° plies. Transverse and longitudinal cracks are assumed to be spaced uniformly and to span the full thickness and width of the 90° and 0° plies, while delaminations are assumed strip-shaped. Spacings between longitudinal and transverse cracks are denoted, respectively, 2*s*_{1} and 2*s*_{2}, while the length of longitudinal and transverse delaminations are denoted 2ℓ_{1} and 2ℓ_{2}, respectively. A set of Cartesian coordinates with the origin in the centre of the laminate is introduced, with *x*_{1}-axis coinciding with the fibre direction in the 90° lamina and *x*_{3}-axis directed through the laminate thickness. The laminate is subjected to general in-plane biaxial tension ( and ) and shear loading ().

In order to analyse *in situ* constrain effect on the stiffness of a particular cracked lamina, the ECM of the damaged laminate is employed [22–24]. In the ECM laminate, all the laminae below and above the damaged lamina under consideration are replaced with homogeneous layers (I and II) having the equivalent constraining effect (figure 2). In-plane stiffness properties of the equivalent constraint layer can be obtained from the laminated plate theory once their stresses and strains are known from micromechanical analysis. Theoretically, the ECM does not impose any restrictions onto the laminate lay-up, and the approach was applied to analysis of quasi-isotropic laminate with matrix cracking in all but 0° layers by Zhang & Herrmann [25].

Application of the ECM approach to cross-ply laminate damaged by transverse and longitudinal matrix cracks and transverse and longitudinal crack-induced delaminations is schematically shown in figure 3. Instead of considering the damaged laminate configuration shown in figure 1, the following two ECMs are analysed instead. In ECM1 (figure 3*a*), the 0° lamina (layer 1) contains damage explicitly, while the 90° lamina (layer 2), damaged by transverse cracking and transverse delaminations, is replaced with the homogeneous layer with reduced stiffness properties. Likewise, in ECM2 (figure 3*b*), the 90° lamina (layer 2) is damaged explicitly, while the damaged 0° lamina is replaced with the homogeneous layer with reduced stiffness. All the quantities associated with the 0° lamina (layer 1) will be henceforth denoted by a sub- or superscript (1), whereas those associated with the 90° lamina (layer 2) with a sub- or superscript (2).

The reduced stiffness properties of the *μ*th layer (*μ*=1,2) damaged by transverse cracking and transverse delaminations (if *μ*=2) or splitting and longitudinal delaminations (if *μ*=1) can be calculated from the laminated plate theory, provided stresses and strains in the explicitly damaged *μ*th layer are known from the analysis of the ECM*μ* (i.e. ECM1 if *μ*=1 and ECM2 if *μ*=2). The reduced elastic properties of the equivalently constraining layer *κ*,*κ*≠*μ*, required in the analysis of the ECM*μ* are supposed to be determined from the analysis of the ECM*κ*. Thus, the problems for ECM1 and ECM2 are inter-related, damage coupling effect is included in the residual stiffness analysis.

## 3. Stress analysis

Owing to the periodicity of damage configuration in the ECM*μ*, only their representative segments (figure 3), containing either a pair of splits or a single transverse crack as well as two strip-shaped delaminations, need to be considered. As the representative segments are symmetric with respect to the mid-plane and their material and geometry are notably uniform in the direction perpendicular to the *x*_{μ}0*x*_{3} plane, the analysis can be further restricted to one quarter of the representative segments. The representative segments of ECM1 and ECM2 can be segregated into perfectly bonded (ℓ_{μ}<|*x*_{μ}|<*s*_{μ}) regions and locally delaminated (|*x*_{μ}|<ℓ_{μ},*μ*=1,2) regions, with no frictional contact between the layers in the latter.

In the perfectly bonded regions (ℓ_{μ}<|*x*_{μ}|<*s*_{μ}) of the ECM*μ*, stresses can be determined from the equilibrium equations
3.1Here are the peak shear stresses at the (0/90) interface of the ECM*μ* in the *x*_{μ}0*x*_{3} plane; are the in-plane microstresses in the *k*th layer of the ECM*μ*, i.e. the stresses averaged across the thickness of the *k*th layer and the width of the ECM*μ* as follows:
3.2In the locally delaminated region (|*x*_{μ}|≤ℓ_{μ}) of the ECM*μ*, the in-plane microstresses in the explicitly damaged *μ*th layer vanish, i.e.
3.3The in-plane microstresses are related to the total stresses applied to the laminate by the following equilibrium equations:
3.4It is assumed that both the explicitly damaged and the equivalently constraining laminae in the ECM*μ* are homogeneous orthotropic, and their constitutive equations, in terms of the in-plane microstresses and microstrains, can be written as
3.5*a*and
3.5*b*where denotes the in-plane stiffness matrix of the explicitly damaged *μ*th layer (a circumflex () is used for representing the elastic properties of the undamaged material) and [*Q*^{(κ)}],*κ*≠*μ*, denotes the in-plane stiffness matrix of the homogeneous orthotropic material of the equivalently constraining *κ*th layer. The in-plane constitutive equations can also be written in terms of strains as
3.6*a*and
3.6*b*where , denote the in-plane compliance matrices of the explicitly damaged *μ*th layer and equivalently constraining *κ*th layer, respectively.

In order to determine the in-plane microstresses in the perfectly bonded region from the equilibrium equations, equation (3.1), the interface shear stresses are expressed in terms of in-plane displacements . Here, it is assumed that the out-of-plane shear stresses , vary linearly with *x*_{3}, which corresponds to a parabolic variation of the in-plane displacements. Besides that, it is assumed that in the 0°-lamina linear variation of the out-of-plane shear stresses is restricted to the region of about one ply thickness (i.e. the nominal thickness of the pre-preg used to make the laminate). We assume that all layers of the laminate have thicknesses in the multiples of the nominal ply thickness. For laminates with thick 0° layer, this appears to offer a more reasonable description of the cracked laminate behaviour. Thus, here the out-of-plane shear stresses are assumed to vary as follows:
3.7where *h*_{s} is the thickness of the shear layer, *m*_{s} is the number of plies in the shear layer and *t* is the ply thickness. After some mathematical calculations and equation rearrangements (e.g. appendix A in [26]), the interface shear stresses are obtained as
3.8where the shear lag parameters *K*_{j} are functions of ply properties
3.9Here, , are the out-of-plane shear moduli of the *k*th layer. As the presence of aligned microcracks does not affect the value of the out-of-plane shear moduli (this fact is emphasized by marking them with a circumflex ()), the shear lag parameters *K*_{j} are the same for ECM1 and ECM2.

The equilibrium equations, equation (3.1), along with expressions for the interface shear stresses (equation (3.8)), the laminate equilibrium equations (equation (3.4)), and constitutive equations (equation (3.6)), provide a full set of equations, which are required for determining the in-plane microstresses in the perfectly bonded regions of the representative segment of the ECM*μ*. For instance, can be found from the following set of eight equations with respect to eight variables:
3.10*a*
3.10*b*
3.10*c*
3.10*d*
3.10*e*
3.10*f*After some rearrangement, this and other similar sets of equations can be reduced to the single differential equations
3.11*a*
3.11*b*where are the laminate constants depending on the layer compliances shear lag parameters *K*_{j} and the layer thickness ratio *χ*=*h*_{1}/*h*_{2}. In detail, they are presented in appendix B of [26]. Given the stress-free boundary conditions at the crack/split surfaces, solutions to equations (3.11) are
3.12*a*and
3.12*b*where *s*_{μ} is crack/split half-spacing and ℓ_{μ} is crack/split tip delamination half-length (figures 1 and 3). Once the in-plane microstresses, equation (3.12), in the explicitly damaged *μ*th layer of the EMC*μ* are known, the laminate macrostresses can be found as
3.13The reduced stiffness properties of the layer *μ*, damaged by transverse cracking or splitting and delaminations, can be determined by applying the laminate plate theory to the EMC*μ* after replacing the explicitly damaged layer with an equivalent homogeneous one. The constitutive equations for the homogeneous layer equivalent to the explicitly damaged *μ*th layer are
3.14Where in order to satisfy compatibility, the macrostrains are assumed to be
3.15

## 4. Stiffness of a damaged lamina

The in-plane reduced stiffness matrix [*Q*^{(μ)}] of the homogeneous layer equivalent to the *μ*th layer of the ECM*μ* is
4.1
4.2

The *in situ* damage effective functions (IDEFs) introduced in [21–23] can be expressed in terms of macrostresses and macrostrains in the *μ*th layer of the ECM*μ* as
4.3*a*and
4.3*b*On substituting macrostresses, calculated from equations (2.13), and macrostrains, calculated from equation (3.15), into equation (4.1), the closed-form expressions for IDEFs are obtained. They represent as functions of relative cracking/splitting density relative delamination area the layer compliances , shear lag parameters *K*_{j} and the layer thickness ratio *χ*:
4.4In detail, the closed-form expressions for the IDEFs for the *μ*th layer of the ECM*μ* are
4.5*a*and
4.5*b*where the constants , *i*=1,2 (appendix C of [26]) depend solely on the layer compliances shear lag parameters *K*_{j} and the layer thickness ratio *χ*. The modified compliances , of the equivalently constraining’ *κ*th layer of the ECM*μ* are determined from the analysis of the ECM*κ* and therefore are functions of the IDEFs Thus, the IDEFs for the *μ*th layer depend implicitly on the damage parameters associated with the *κ*th layer.

The IDEFs for both layers form a system of simultaneous nonlinear algebraic equations
4.6*a*and
4.6*b*This system is solved computationally using a direct iterative procedure. It is carried out in such a way that the newly calculated IDEFs are used to evaluate the reduced stiffnesses of the equivalently constraining *κ*th layer repeatedly until the difference between two iterative steps meets the prescribed accuracy. Consequently, all four IDEFs , *k*=1,2 are determined as functions of damage parameters If interactions between damage modes in different laminae are neglected, IDEFs associated with the *μ*th layer will depend only on damaged parameters for that layer.

Verification of the ECM/2-D shear lag approach in absence of crack-induced delaminations was carried out in [15,17,19,27]. After comparison with other existing models by Hashin [28], Tsai & Daniel [29] and Henaff-Gardin *et al.* [30,31] describing stiffness reduction of carbon-fibre-reinforced polymer and glass-fibre-reinforced polymer cross-ply laminates due to transverse cracking and splitting, the following conclusions were reached in [15]. As far as the reduction of Young’s modulus is concerned, the ECM/2-D shear lag approach is in very good agreement with other models. Its predictions are closer to the lower bound established by Hashin [28] than the results of Henaff-Gardin *et al.* [30]. For Poisson’s ratio, the ECM/2-D shear lag approach predictions are close to those of Henaff-Gardin *et al.* [30], although for small values of the damage parameter (relative crack/split spacing) the reduction predicted by the ECM/2-D shear lag approach is greater than of Henaff-Gardin *et al.* [30]. Predictions based on the variational approach of Hashin [28] are far away from these results. The shear modulus reduction ratio predicted by Tsai & Daniel [29] is, in most cases, within 10% of the ECM/2-D shear lag approach value.

It is worth mentioning here that the model of Tsai & Daniel [29] and the present ECM/2-D shear lag approach yield exactly the same analytical expression for the shear modulus reduction ratio due to transverse cracking, if the thickness of the shear layer in the ECM/2-D shear lag approach is taken equal to that of the 0^{o} lamina, i.e. if *h*_{s}=*h*_{1}:
4.7For transverse cracking combined with splitting, Tsai & Daniel [29] suggested a semi-empirical expression for the shear modulus reduction ratio based on the ‘superposition’ of solutions for a single set of cracks as
4.8The value of the shear modulus reduction ratio obtained by Tsai & Daniel [29] using the finite difference iteration appeared to be within 1% of the value given by equation (4.8). The present ECM/2-D shear lag model, if the interaction between transverse cracks and splits is neglected and the shear layer has the thickness of the 0° lamina, yields an expression
4.9It may be seen from equations (4.8) and (4.9) that the two expressions differ by the underlined terms and *ρ***G*≤*ρ*_{G}. In the absence of splitting ) they are both reduced to equation (4.7). In some cases, though, the error of the semi-empirical expression, equation (4.8), suggested by Tsai & Daniel [29] can be as big as 20%. The ECM/2-D shear lag approach is in good agreement with the results of Henaff-Gardin *et al.* [30] for the shear modulus reduction.

## 5. Results and discussion

Stiffness degradation in cross-ply laminates due to different damage modes and their combinations is examined below. All results given below were obtained taking into account the interaction between damage modes in the adjacent layers. Up to 12 iterations are required to solve a set of simultaneous nonlinear equations, equations (4.6), with an accuracy of 10^{−9}. The number of iterations increases along with the crack density and relative delamination area.

Figure 4 shows stiffness degradation in E-glass/LY556 epoxy [26] [0/90/0] and [0/90_{8}/0] cross-ply laminates as a function of the transverse crack density *C*_{2} in the 90° layer. The layer thicknesses *h*_{1} and *h*_{2} are determined from the laminate lay-up; thickness of the shear layer is taken as *h*_{s}=*t*. Longitudinal Young’ modulus, shear modulus and major Poisson’s ratio are normalized by their values in the undamaged state. As can be seen from figure 4*a*,*b*, all these properties undergo degradation as the matrix crack density increases, with Poisson’s ratio appearing to be the most affected by transverse cracking. The thickness of the 90° layers plays an important role, since the thicker the 90° layer, the bigger is the reduction observed. Transverse ply thickness and the thickness ratio of 90° layer to constraining 0° layers are the important parameters controlling resistance to matrix cracking. Zhang *et al*. [23] proposed to use a resistance curve, analogous to the *R*-curve concept of classical fracture mechanics, as a measure of the composite resistance to crack initiation and growth
5.1where *G* is the strain energy release rate associated with matrix cracking, *G*_{R} is the laminate resistance to matrix cracking, *G*_{IC} is the critical energy release rate for damage nucleation, and *G*_{0} and *R* are laminate constants. Parameters *G*_{IC}, *G*_{0} and *R* are not independent of stacking sequence, but remain constant as long the thickness ratio of the constraining layer to 90° remains the same.

When a cross-ply laminate is subjected to biaxial loading, matrix cracking may occur concurrently in both plies leading to formation of transverse and longitudinal matrix cracks. The combined effect of these cracks on stiffness properties of [0/90/0] laminate is shown in figure 5 for the case when the longitudinal and transverse crack densities are equal.

In cross-ply laminates with thick 90° layer subjected to uniaxial loading, strip-shaped delaminations begin to initiate and grow from the tips of matrix cracks at the 0°/90° interface. The effect of these delaminations on stiffness properties of [0/90_{8}/0] laminate is shown in figure 6 as a function of relative delamination area. Transverse crack density is taken as 2 cracks cm^{−1}, and the values of normalized stiffness properties for *D*_{2}=0 correspond to stiffness degradation due to matrix cracking without delamination. It can be seen from figure 6 that crack tip delamination contributes significantly to stiffness degradation of the laminate, and therefore has to be taken into account in the post-initial failure models.

Figure 7 shows stiffness degradation in T800H/3631 carbon/epoxy [0/90_{n}]_{s}, *n*=2,4,6, cross-ply laminates containing transverse cracks and delaminations. Longitudinal Young’s modulus, in-plane shear modulus and Poisson’s ratio, normalized by their values in the undamaged state, are plotted as a function of transverse crack density. The relative delamination area is 10%, which corresponds to ℓ/*s*=0.1. For the axial modulus, predictions are compared to experimental data obtained by Takeda & Ogihara [32] and appear to be in good agreement. However, predictions show that reduction in shear modulus and Poisson’s ratio due to crack tip delamination is more significant.

Henaff-Gardin *et al.* [31] observed damage development in cross-ply [0_{4}/90_{4}]_{s} T300/914 carbon/epoxy laminates during thermal cycling. The cycle consisted of cooling to −200°C and heating to +90°C. Crack density in the 90° and 0° plies was measured; however the size of growing delaminations that accompanied longitudinal cracks was not. In figure 8*a*, predictions of reduction in the longitudinal Young’s modulus, shear modulus and Poisson’s ratio, normalized by their values in the undamaged state, are shown along with the measured crack densities in the 90° and 0° plies as a function of number of cycles. As cracks develop, the shear modulus and Poisson’s ratio undergo significant reduction, while reduction in axial modulus remains less than 5%. This indicates that the shear modulus and Poisson’s ratio could be much better parameters to characterize stiffness degradation of the laminate than the longitudinal modulus.

Since the size of the delamination area was not measured during cycling, reduction of stiffness properties of [0_{4}/90_{4}]_{s} T300/914 laminate due to delaminations was predicted using assumed delamination sizes. Strip-width of the transverse delamination was set to zero, while that of the longitudinal delamination was allowed to vary from zero to 50%. In other words, longitudinal delaminations were assumed to have propagated from the crack tip to one quarter of the distance between two cracks. This seems to be a reasonable assumption, consistent with X-ray radiographs obtained by Henaff-Gardin *et al.* [31]. In figure 8*b*,*c*, predicted reductions of the longitudinal, transverse and shear moduli as well as Poisson’s ratio, normalized by their values in the undamaged state, are plotted as a function of the relative delamination area. The axial modulus appears to be unaffected by the growth of delamination, while transverse modulus is further reduced, but not significantly (figure 8*b*). The reduction in the shear modulus is more pronounced than in Poisson’s ratio (figure 8*c*). Crack densities in 90° and 0° plies were taken as *C*_{2}=4.5 cracks cm^{−1} and *C*_{1}=3 cracks cm^{−1}, respectively, which corresponds to saturation values reached during −200°C/+90°C cycling. Under uniaxial loading, longitudinal delaminations appear to be more important than the transverse ones, since they result in isolation of the portions of the load-bearing 0° plies, which become prone to fibre breakage. Under biaxial loading, the importance of one set of delaminations over the other depends very much on the biaxiality and ply thickness ratios.

## 6. Conclusion

Although the approach described in this paper has not attempted to predict ultimate laminate failure, it does present a methodology for predicting degraded stiffness properties of the laminae and hence the laminate, in the case when there are various kinds of intra- and interlaminar damage interacting with each other in the same and/or adjacent plies of the laminate. The approach is based on the ECM of the damaged laminate and takes into account damage mode interaction. Our predictions show that the effect of longitudinal matrix cracking is more pronounced on Poisson’s ratio than on the shear modulus; however, the reduction in the shear modulus due to transverse delamination is the most significant when compared with the reduction observed in the axial or transverse elastic moduli.

Theoretically, ECM does not impose any restrictions on the laminate lay-up, and the approach based on ECM was successfully applied to the prediction of degraded stiffness properties due to matrix cracking in all but 0° layers of quasi-isotropic laminates. It should be noted that for the model to be applied the type, location and amount of damage present need to be specified. For this accurate and reliable structural health monitoring techniques are urgently required [33–35]. Also the triggering of resin cracking and delamination could be delayed to higher applied loads if tougher resin systems are employed [36].

## Authors' contributions

M.K. and C.S. performed modelling, carried out analysis of results and drafted the paper. Both authors gave final approval for publication.

## Competing interests

We declare we have no competing interests.

## Funding

Financial support of this research by the Engineering and Physical Sciences Research Council (EPSRC/GR/L51348) and the British Ministry of Defence is gratefully acknowledged.

## Footnotes

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

- Accepted February 22, 2016.

- © 2016 The Author(s)

Published by the Royal Society. All rights reserved.