Damage history in composite laminates: Matrix cracks leading to delaminations

Abstract

Strain energy release rate calculations for various cases of delaminations emanating from matrix cracks are developed and used to predict the onset of delaminations and their growth size as a function of applied tension and shear loads in composite laminates. The method determines the matrix crack spacing, the delamination onset load, the delamination size at onset and, through the use of a newly proposed delamination resistance curve, the size of delaminations as they grow under load. The method can be applied to any symmetric laminate. Comparisons to test results in the literature for a variety of layups and materials shows very good agreement with the exception of cases where significant edge delaminations appear before delaminations caused by matrix cracks.

Introduction

In broad terms, there are three major types of damage in composites: Broken fibers, matrix cracks and delaminations. The discussion in this work focuses on the interaction of delaminations with matrix cracks. It is well-known^1–3 that, often, delaminations start at the tips of matrix cracks parallel to the fibers. In a ply with matrix cracks, the load is transferred to adjacent plies around these cracks and interlaminar stresses develop. These may cause delaminations at the interface between the cracked ply and its neighbours. The presence of matrix cracks and delaminations reduces the stiffness properties of the laminate and may contribute to pre-mature failure.^4,5

Various methods to analytically predict the creation of matrix cracks have been proposed over the years: Continuum mechanics methods, [for example Ref. 6,7], shear lag analyses, [for example Ref. 8,9], and variational approaches.^10,11 Talreja ⁶ derived constitutive relations describing the behaviour of a laminate with different damage modes present, which can be used to describe damage evolution under static or fatigue loads. Varna et al ⁷ focused on laminates containing 90° plies and proposed an approach to predict matrix cracks and their effect on stiffness of such laminates. Kashtalyan and Soutis, ⁸ combined matrix cracks in angle-ply laminates with delaminations. Nairn, ⁹ examined the use of shear lag methods for analysing unidirectional laminates. Hashin,^10,11 was among the first to use a variational approach to determine stresses in cracked laminates and their effect on stiffness of cross-ply laminates. In addition to the above, Socci and Kassapoglou ⁵ solved the governing equations to determine matrix crack creation and multiplication for laminates under tension and shear loads.

Delaminations in laminates under in-plane loads, typically form at free edges because of mismatch in stiffness properties of the plies^12,13 or at tips of matrix cracks. Of particular interest here are the latter where the interlaminar stresses that develop at the tips of matrix cracks may exceed the strength of the material in the resin layer between plies and lead to delaminations. Nairn and Hu ¹⁴ used a variational approach to obtain expressions for the strain energy release rate expressions to predict the critical crack spacing for delamination onset for a class of laminates containing 90° plies. A similar approach was proposed by Lem and Lui ¹⁵ where they compared the energy released when cracks or delaminations form. A three-dimensional finite element analysis was used by Salpekar and O’Brien ¹⁶ to study the formation of triangular delaminations at matrix crack tips. Takeda and Ogihara, ³ obtained expressions for the onset and growth of delaminations at matrix crack tips and compared to test results. Finally, Carraro et al, ² developed expressions for the stresses around the matrix crack tips using generalized stress intensity factors and proposed a fracture criterion using a critical value of the mode I Generalized Stress Intensity Factor.

The complexity of the stress field at matrix crack tips and around delaminations emanating from these crack tips makes the development of analytical methods difficult and, often, limited to a class of laminates and/or applied loading situations. Also, localized defects such as voids or matrix rich regions accelerate or delay the creation of cracks and delamination. Furthermore, measuring the crack density exactly poses problems because of the difficulty of pin-pointing the exact applied load at which new cracks form. Finally, as has been observed by several researchers,^17,18 the symmetry of delamination creation assumed in many models breaks down sometimes as some delaminations forming in one ply interface stop growing and/or jump to other ply interfaces.

The present communication builds on previous work and proposes a fracture mechanics-based approach to analytically predict the onset and growth of delaminations emanating from matrix crack tips for symmetric composite laminates under tension or shear loads.

Approach

Matrix crack formation

The starting point is a laminate in some of the plies of which matrix cracks developed under some type of loading, see Figure 1. The layup can be arbitrary as long as it is symmetric, and the loading can be any combination of in-plane tension and shear. Subsequent loading may either increase the density of these cracks or/and lead to the creation of delaminations. In order to model the delamination onset at tips of matrix cracks, it is critical that the crack density in the cracked plies be correctly predicted. Otherwise, the stiffness changes caused by the presence of these cracks is not correctly accounted for. Here, the approach developed by Socci and Kassapoglou ⁵ is used. For convenience, the important equations from that work are briefly repeated here.

Figure 1.

Symmetric laminate with ply with matrix cracks under in-plane loading.

Given a ply with matrix cracks with spacing D as shown in Figure 1, the reduced transverse Young’s modulus (perpendicular to the fibers) in the cracked ply is given by: ⁵

$\frac{E_{2red}}{E_2}=1+\frac{8}{{\pi }^2}\sum _n\frac{1}{n^2}\frac{2e^{\frac{{\phi }_{n}D}{2}}-2e^{-\frac{{\phi }_{n}D}{2}}+{\phi }_{n}D\left(e^{\frac{{\phi }_{n}D}{2}}+e^{-\frac{{\phi }_{n}D}{2}}\right)}{e^{-{\phi }_{n}D}-e^{{\phi }_{n}D}-2{\phi }_{n}D}$

(1)where E₂ is the transverse Young’s modulus for the uncracked ply, φ_n = nπ/t_c with t_c the thickness of the cracked ply and D the crack spacing. By Maxwell’s relation, the minor Poisson’s ratio of the cracked ply is given by:

${\nu }_{21red}=\frac{E_{2red}}{E_1}{\nu }_{12}$

(2)with E₁ and ν₁₂ the Young’s modulus parallel to fibers and the major Poisson’s ratio of the pristine ply.

If the ply is, locally, under shear and, accounting for non-linearity of the shear stress-strain curve of the material, the shear modulus of a cracked ply is given by: ⁵

$\frac{G_{12red}}{G_{12}}=\frac{1-8B_n}{1-\frac{8B_n\left(1-\frac{{\gamma }_y}{{\gamma }_a}\right)\left[\left(1-8B_n\right)\left(1-\frac{{\gamma }_y}{{\gamma }_a}\right)+\frac{{\gamma }_y}{{\gamma }_a}\right]}{1-8B_n+\frac{{\gamma }_y}{{\gamma }_a}}}$

(3)where

$B_n=\sum _n\frac{e^{k_n\frac{D}{2}}+{e^{k_{n}D}e}^{-k_n\frac{D}{2}}}{{\left(n\pi \right)}^2\left(1+e^{k_{n}D}\right)}$

(4)with γ _y the yield shear strain of the material, k _n = nπ/t _c and n odd. It is interesting to note that, unlike the reduced transverse modulus E _2red , the reduced shear modulus is a function of the local applied strain magnitude γ _a .

The load at which matrix cracks will multiply is found using a criterion based on the average energy density between cracks. ⁵ For an arbitrary combination of transverse and shear strains applied locally to the ply, a new crack will appear half-way between existing cracks when:

$\frac{\Delta {\sigma }^2}{2E_2}+\frac{\Delta {\tau }^2}{2G_{12}}=\Delta U_{davcrit}$

(5)where Δσ and Δτ refer to the normal and shear stresses in the ply before and after a new crack is created and ΔU _davcrtit is the thickness averaged critical energy density for the given combination of normal and shear strains. ⁵ This critical energy density changes with combination of σ and τ and is determined by considering the limiting case of large crack spacing D where there is enough distance between cracks for the ply to develop the pristine values of σ and τ. In such a case, local failure of the ply would coincide with creation of a new crack. A Hashin-type failure criterion is used. Then, ⁵ :

$\Delta U_{davcrit}=\frac{{\sigma }_{yav}^2}{2E_2}+\frac{{\tau }_{xyav}^2}{2G_{12}}$

(6)provided the average transverse stress σ _yav and the average shear stress τ _xyav in the ply satisfy the criterion:

$\frac{{\sigma }_{yav}^2}{{\left(Y_{is}^t\right)}^2}+\frac{{\tau }_{xyav}^2}{{\left(S_{is}\right)}^2}=1$

(7)

As the above equations will form the starting point for the determination of delamination onset and growth it is important to establish confidence in them. This is done by comparing to published experimental results. The size of the crack spacing in 90° plies as a function of applied laminate stress for [0_n/90_m]s and [±25/90_n]s laminates for different values of n and m and for different materials is compared to test results from^2,19–22 in Figure 2.

Figure 2.

Comparison of predicted crack spacing to tests for various layups and materials. (a) UT-E500/Epikote RIMR235 fiberglass with tests from Ref. 2. (b) AS/3501-6 Carbon/Epoxy with tests from Refs.20,21. (c) AS4/3501-6 Carbon/Epoxy with tests from Refs.22,23. (d) IM6/Avimid K Carbon/Epoxy with tests from Refs.22,23.

As is seen from Figure 2, the predictions are in good to excellent agreement with the test results. One interesting observation is that for the [0/90₂]s laminate in Figure 4(a) and for the four laminates in Figure 4(c) and (d), the curves and test data seem to asymptote to an horizontal line. This corresponds to crack saturation where further load increase does not cause the creation of any more matrix cracks. This would correspond to the point where delaminations may start at the tips of these cracks. The precise crack spacing at which this happens will be determined in the next section.

Delamination onset – delaminations emanating from matrix cracks

 The results presented in the previous section establish enough confidence in the accuracy of the method in predicting matrix crack spacing as a function of applied load. This forms the starting point for determining when delaminations will start at the matrix crack tips. First, the strain energy release rate corresponding to a given distribution of matrix cracks is determined. The strain energy of the laminate under tension and shear shown in Figure 1 can be written as:

$U=\frac{1}{2}\iiint \left[E_x^c{{\epsilon }_x}^2+{g_{xy}^c{\gamma }_{xy}}^2\right]dV=\frac{1}{2}\left[E_x^c{{\epsilon }_x}^2+{g_{xy}^c{\gamma }_{xy}}^2\right]Lwh$

(8)where ε _x and γ _xy are laminate axial and shear strains, E _x ^c and g _xy ^c are laminate Young’s modulus and shear modulus accounting for the presence of matrix cracks and L, w, and h are the length, width and height of the laminate respectively. Note that “g” is used for shear modulus instead of “G” to avoid confusion with the symbol for the strain energy release rate. The strain energy release rate relates the rate at which the energy is dissipated by a material when a crack forms. In this case, the right hand side of equation (8) must be differentiated with respect to the crack area which, for a cracked ply (or multiple plies with same orientation) of (total) thickness t_c is:

$A_c=t_{c}w\frac{L}{D}$

(9)

Then, the strain energy release rate for matrix cracking is given by:

$G^c=\frac{\partial U}{\partial A_c}=\frac{1}{2}\frac{h}{t_c}\left[\frac{\partial E_x^c}{\partial \left(\frac{1}{D}\right)}{\epsilon }_x^2+\frac{\partial g_{xy}^c}{\partial \left(\frac{1}{D}\right)}{\gamma }_{xy}^2\right]$

(10)

In order to proceed, the dependence of E _x ^c and g _xy ^c on 1/D (matrix crack density) must be determined. This follows the standard approach for determining the engineering constants of a symmetric laminate using the classical laminated-plate theory. In local coordinates, parallel and perpendicular to the fibers, the stiffness values for a cracked ply with matrix cracks are given by

$Q_{11}=\frac{E_1}{\left(1-{\nu }_{12}{\nu }_{21red}\right)}=\frac{E_1}{\left(1-{\nu }_{12}{\nu }_{21}\frac{E_{2red}}{E_2}\right)}$

(11)

$Q_{22}=\frac{E_{2red}}{\left(1-{\nu }_{12}{\nu }_{21red}\right)}=\frac{E_{2red}}{\left(1-{\nu }_{12}{\nu }_{21}\frac{E_{2red}}{E_2}\right)}$

(12)

$Q_{12}=\frac{{\nu }_{12}{E}_{2red}}{\left(1-{\nu }_{12}{\nu }_{21red}\right)}=\frac{{\nu }_{12}{E}_{2red}}{\left(1-{\nu }_{12}{\nu }_{21}\frac{E_{2red}}{E_2}\right)}$

(13)

$Q_{66}=G_{12red}$

(14)where equation (2) was used to obtain the final expression in the right hand sides and E _2red and G _12red are given by equations (1) and (3). Denoting then, d(..)/d(1/D) by (…)’ the derivatives needed to obtain the strain energy release rate are obtained:

${Q_{11}}^{\prime }=\frac{{{\nu }_{12}}^2{E}_{2red}^{\prime }}{{\left(1-{\nu }_{12}{\nu }_{21}\frac{E_{2red}}{E_2}\right)}^2}$

(15)

${Q_{22}}^{\prime }=\frac{E_{2red}^{\prime }}{1-{\nu }_{12}{\nu }_{21}\frac{E_{2red}}{E_2}}\left[1+\frac{{{\nu }_{12}}^2\frac{E_{2red}}{E_1}}{1-{\nu }_{12}{\nu }_{21}\frac{E_{2red}}{E_2}}\right]$

(16)

${Q_{12}}^{\prime }=\frac{{{\nu }_{12}E}_{2red}^{\prime }}{1-{\nu }_{12}{\nu }_{21}\frac{E_{2red}}{E_2}}\left[1+\frac{{{\nu }_{12}}^2\frac{E_{2red}}{E_1}}{1-{\nu }_{12}{\nu }_{21}\frac{E_{2red}}{E_2}}\right]$

(17)

${Q_{66}}^{\prime }=G_{12red}^{\prime }$

(18)

In addition,

$E_{2red}^{\prime }=\frac{\partial E_{2red}}{\partial \left(1/D\right)}=-\frac{4D^2{E}_2}{t_c^2}\sum _n\frac{e^{\frac{{\phi }_{n}D}{2}}\left(e^{{\phi }_{n}D}-1\right)\left[\frac{1}{D}\left(e^{{2\phi }_{n}D}-1\right)-2{\phi }_n{e}^{{\phi }_{n}D}\right]}{{\left[2{\phi }_n{e}^{{\phi }_{n}D}+\frac{1}{D}\left(e^{{2\phi }_{n}D}-1\right)\right]}^2}$

(19)

$G_{12red}^{\prime }=\frac{\partial G_{12red}}{\partial \left(1/D\right)}=\frac{-8B_n^{\prime }}{1-\frac{8B_n\left(1-\frac{{\gamma }_y}{{\gamma }_a}\right)\left[\left(1-8B_n\right)\left(1-\frac{{\gamma }_y}{{\gamma }_a}\right)+\frac{{\gamma }_y}{{\gamma }_a}\right]}{1-8B_n+\frac{{\gamma }_y}{{\gamma }_a}}}+\frac{64\left(1-B_n\right)B_n\left(1-\frac{{\gamma }_y}{{\gamma }_a}\right)\left[\left(1-8B_n\right)\left(1-\frac{{\gamma }_y}{{\gamma }_a}\right)+\frac{{\gamma }_y}{{\gamma }_a}\right]B_n^{\prime }}{{\left[1-\frac{8B_n\left(1-\frac{{\gamma }_y}{{\gamma }_a}\right)\left[\left(1-8B_n\right)\left(1-\frac{{\gamma }_y}{{\gamma }_a}\right)+\frac{{\gamma }_y}{{\gamma }_a}\right]}{1-8B_n+\frac{{\gamma }_y}{{\gamma }_a}}\right]}^2{\left(1-8B_n+\frac{{\gamma }_y}{{\gamma }_a}\right)}^2}$

(20)with

$B_n^{\prime }=\sum _n\frac{{k_{n}e}^{k_n\frac{D}{2}}\left(e^{k_{n}D}-1\right)}{{\left(n\pi \right)}^2{\left(1+e^{k_{n}D}\right)}^2}{D}^2$

(21)

For a symmetric laminate, E _x ^c and g _xy ^c can be calculated using the A matrix. The results can then be differentiated with respect to 1/D with the use of equations (19) and (20) to get the strain energy release rate G^c for the laminate with cracked plies. An example of how the value of G ^c changes with crack spacing for different thickness of cracked plies in various laminates under tension is shown in Figure 3.

Figure 3.

Energy release rate dependence on matrix crack spacing for various layups. Applied strain = 0.05.

As the applied load increases from zero, the (average) crack spacing, which is large and non-uniform at the beginning, will decrease which means one moves from right to left in Figure 3. For low applied load and, hence, large values of D, the creation of an additional crack is not affected by existing cracks. There is no interaction and all cracks have the same size. The change in strain energy for a given constant crack size will be the same every time a new crack forms. As a result, the strain energy release rate will be essentially constant for large values of D. As the applied load increases and the crack spacing D decreases, a point is reached where crack spacing is uniform and a new crack would form at the mid-point between cracks. ⁵ This was shown in reference 5 to correspond to D = 2.966 t _c . The amount of energy released when a crack forms increases and the energy release rate goes through a maximum. It should be noted that, other than the case θ = 90 ^o for the [0/θ₂]s layups in Figure 3, all laminates have both tension and shear strains present in the ply with θ orientation. As the value of θ decreases from 90^o the magnitude of the shear stress in the θ ply increases and that of the transverse tension stress decreases until θ = 45^o. As was found in Ref. 5, a higher shear strain magnitude is needed to initiate a new crack than a tension strain. Consequently, as can be seen from Figure 3, laminates with shear stress dominant in the θ ply, are less prone to matrix crack creation and thus the strain energy release rate is lower. For θ < 45°, both the magnitudes of shear and transverse stress in the θ ply decrease leading to further reduction in the value of G ^c .

Using equation (10), the effect of various parameters on the strain energy release rate was found. The basic trends are as follows: (a) G ^c increases significantly when t _c , the thickness of the cracked ply, increases, (b) Increasing t _c shifts the G ^c versus D graph to the right, towards larger crack spacings, (c) Increasing E ₁ /E ₂ increases G ^c , (d) G^c increases by a small amount when t ₁ , the thickness of the ply without cracks increases. (d) The point at which G ^c departs significantly from being constant depends only on t _c . At this point, the onset of delamination from existing matrix cracks can be examined. A variety of possibilities, depending on the resulting pattern is examined and is shown in Figure 4. In all cases, the delaminations shown in blue color extend across the laminate width (perpendicular to the page in Figure 4). It is understood that in a real specimen, combinations of the above may appear as a result of local defects such as waviness, porosity and non-uniform matrix content. The cases selected here serve as an indicator of what to expect in reality and will provide a range of results. The approach is analogous to that used by O’Brien ¹³ for edge delaminations and makes use of equations (8) and (10) but now the effect of the delaminations on the laminate stiffness must be calculated. Focusing on Case 1 of Figure 4, the laminate can be divided in repeating units as shown in Figure 5.

Figure 4.

Possible configurations of delaminations emanating from matrix cracks.

Figure 5.

Repeating unit to obtain longitudinal stiffness of a laminate with cracks and delaminations (case I).

For simplicity in showing the equations and understanding how different variables contribute, the case of a cross-ply laminate will be assumed here but the results will be generalized to any symmetric laminate under tension and shear. Referring to Figure 5:

$E_{in}=\frac{1}{1-{\nu }_{12}{\nu }_{21red}}\left[E_{2red}-{\nu }_{12}^2\frac{E_{2red}^2}{E_1}\right]$

(22)

$E_{out}=\frac{1}{1-{\nu }_{12}{\nu }_{21}}\left[E_1-{\nu }_{12}^2{E}_2\right]$

(23)

Strain compatibility of the three sub-laminates with the delaminations present leads to:

$E_x^{del}=\frac{2t_1{E}_{out}+t_c{E}_{in}}{2t_1+t_c}$

(24)

The stiffness of the three sub-laminates on the right of Figure 5 is that of a laminate with matrix cracks only, E _x ^c which was calculated before. Recognizing that the entire laminate consists of the pattern of Figure 5 repeating along the length of the laminate and that each repeating unit has the same applied axial force per unit width N _x , the overall axial stiffness of the laminate is obtained treating the repeating units as springs in series:

$E_x^{c+d}=\frac{D{E}_x^{del}{E}_x^c}{a{E}_x^c+\left(D-a\right)E_x^{del}}$

(25)where a is the delamination length (parallel to the loading shown in Figure 4) and the superscript “c + d” denotes cracks and delaminations present. In an analogous way but considering shear strains and load per unit width, the shear modulus of the laminate is found:

$g_{xy}^{c+d}=\frac{D{g}_{xy}^{del}{g}_{xy}^c}{a{g}_{xy}^c+\left(D-a\right)g_{xy}^{del}}$

(26)

Case II can be shown to give exactly the same expressions for E_x^c+d and g_xy^c+d. For case III, the pattern to be analysed is shown in Figure 6. An analogous approach leads to the expressions:

$E_x^{c+d}=\frac{2D{E}_x^{del}{E}_x^c}{a{E}_x^c+\left(2D-a\right)E_x^{del}}$

(27)

$g_{xy}^{c+d}=\frac{2D{g}_{xy}^{del}{g}_{xy}^c}{a{g}_{xy}^c+\left(2D-a\right)g_{xy}^{del}}$

(28)

Figure 6.

Repeating unit to obtain longitudinal stiffness of a laminate with cracks and delaminations (case III).

Finally, for case IV, the repeating unit is shown in Figure 7 and the resulting stiffness expressions are:

$E_x^{c+d}=\frac{D{E}_x^{del*}{E}_x^c}{a{E}_x^c+\left(D-a\right)E_x^{del*}}$

(29)

$g_{xy}^{c+d}=\frac{D{g}_{xy}^{del*}{g}_{xy}^c}{a{g}_{xy}^c+\left(D-a\right)g_{xy}^{del*}}$

(30)

Figure 7.

Repeating unit to obtain longitudinal stiffness of a laminate with cracks and delaminations (case IV).

Note that the quantities with superscript “del*” refer to the right most laminate in Figure 7.

The basic stiffness values in equations (25)–(30) need only be calculated once for each stacking sequence as they will not change across cases. The approach just presented assumed only one cracked ply or block of contiguous plies. If failure analysis suggests that more than one plies will develop cracks, the same approach can be used provided that the repeating units and associated derivation include all cracked plies and resulting delaminations.

The strain energy release rate can now be calculated using an equation analogous to equation (10) with appropriate substitution for the stiffness values depending on the case. This energy release rate will describe energy released when the delaminations grow. It is assumed that, once the delaminations appear, further load increase will lead to delamination growth and not the creation of matrix cracks whose density has reached saturation. As a result, the crack area created is the area of the delaminations. This means that equation (9) now changes according to the case to:

$\begin{array}{l}{A}_d=2aw, Case\hspace{0.17em}I\\ A_d=2\left(2\frac{a}{2}\right)w=2aw, Case\hspace{0.17em}II\\ A_d=2aw, Case\hspace{0.17em}III\\ A_d=2\frac{a}{2}w=aw, Case\hspace{0.17em}IV\end{array}$

(31a-d)

Differentiating the strain energy with respect to A_d leads to the following expressions for the strain energy release rate for each of the four cases:

$\begin{array}{l}{G}^{c+d}=\frac{1}{4}h\left(\frac{L}{D}\right)\left[\frac{E_x^c\left(\frac{E_x^c}{E_x^d}-1\right)}{{\left(1-\frac{a}{D}+\frac{a}{D}\frac{E_x^c}{E_x^d}\right)}^2}{\epsilon }_x^2+\frac{g_{xy}^c\left(\frac{g_{xy}^c}{g_{xy}^d}-1\right)}{{\left(1-\frac{a}{D}+\frac{a}{D}\frac{g_{xy}^c}{g_{xy}^d}\right)}^2}{\gamma }_{xy}^2\right],Case\hspace{0.17em}I\\ G^{c+d}=\frac{1}{4}h\left(\frac{L}{D}\right)\left[\frac{E_x^c\left(\frac{E_x^c}{E_x^d}-1\right)}{{\left(1-\frac{a}{D}+\frac{a}{D}\frac{E_x^c}{E_x^d}\right)}^2}{\epsilon }_x^2+\frac{g_{xy}^c\left(\frac{g_{xy}^c}{g_{xy}^d}-1\right)}{{\left(1-\frac{a}{D}+\frac{a}{D}\frac{g_{xy}^c}{g_{xy}^d}\right)}^2}{\gamma }_{xy}^2\right],Case\hspace{0.17em}II\\ G^{c+d}=\frac{1}{2}h\left(\frac{L}{D}\right)\left[\frac{E_x^c\left(\frac{E_x^c}{E_x^d}-1\right)}{{\left(2-\frac{a}{D}+\frac{a}{D}\frac{E_x^c}{E_x^d}\right)}^2}{\epsilon }_x^2+\frac{g_{xy}^c\left(\frac{g_{xy}^c}{g_{xy}^d}-1\right)}{{\left(2-\frac{a}{D}+\frac{a}{D}\frac{g_{xy}^c}{g_{xy}^d}\right)}^2}{\gamma }_{xy}^2\right],Case\hspace{0.17em}III\\ G^{c+d}=\frac{1}{2}h\left(\frac{L}{D}\right)\left[\frac{E_x^c\left(\frac{E_x^c}{E_x^{del*}}-1\right)}{{\left(1-\frac{a}{D}+\frac{a}{D}\frac{E_x^c}{E_x^{del*}}\right)}^2}{\epsilon }_x^2+\frac{g_{xy}^c\left(\frac{g_{xy}^c}{g_{xy}^{del*}}-1\right)}{{\left(1-\frac{a}{D}+\frac{a}{D}\frac{g_{xy}^c}{g_{xy}^{del*}}\right)}^2}{\gamma }_{xy}^2\right],Case\hspace{0.17em}IV\end{array}$

(32a-d)

It can be seen from equations (31a-d) that cases I and II give identical results. Case III has a bigger denominator in the brackets indicating that, for the same values of a and D, the strain energy release rate will be smaller. Case IV is hard to anticipate its magnitude by inspection because of the quantities E _x ^del* and g _xy ^del* whose effect is hard to evaluate without explicit calculations. An approximation to the onset of delamination can be obtained by taking the limit of equations (32a-d) as the delamination size a tends to zero. The resulting expressions are:

$\begin{array}{l}{G}_{onset}^{c+d}=\frac{1}{4}h\left(\frac{L}{D}\right)\left[E_x^c\left(\frac{E_x^c}{E_x^d}-1\right){\epsilon }_x^2+g_{xy}^c\left(\frac{g_{xy}^c}{g_{xy}^d}-1\right){\gamma }_{xy}^2\right],Case\hspace{0.17em}I\\ G_{onset}^{c+d}=\frac{1}{4}h\left(\frac{L}{D}\right)\left[{E_x^c\left(\frac{E_x^c}{E_x^d}-1\right)\epsilon }_x^2+{g_{xy}^c\left(\frac{g_{xy}^c}{g_{xy}^d}-1\right)\gamma }_{xy}^2\right],Case\hspace{0.17em}II\\ G_{onset}^{c+d}=\frac{1}{8}h\left(\frac{L}{D}\right)\left[E_x^c\left(\frac{E_x^c}{E_x^d}-1\right){\epsilon }_x^2+g_{xy}^c\left(\frac{g_{xy}^c}{g_{xy}^d}-1\right){\gamma }_{xy}^2\right],Case\hspace{0.17em}III\\ G_{onset}^{c+d}=\frac{1}{2}h\left(\frac{L}{D}\right)\left[E_x^c\left(\frac{E_x^c}{E_x^d}-1\right){\epsilon }_x^2+{g_{xy}^c\left(\frac{g_{xy}^c}{g_{xy}^{del*}}-1\right)\gamma }_{xy}^2\right],Case\hspace{0.17em}IV\end{array}$

(33a-d)

A plot of equations (33a-d) as a function of crack spacing for a [0/90₂]s UT-E500/Epikote RIMR235 is shown in Figure 8. In addition to the four curves from equations (33a-d), the onset of edge delamination predicted by O’Brien’s method, ¹³ and given by equation (34) below, is also shown as a horizontal line since it is derived for nominally pristine laminate without matrix cracks. In equation (34), E _LAM is the axial stiffness of the pristine laminate and E ^* is the sum of the stiffnesses of the three sub-laminates created by edge delaminations at the 0/90 interfaces. E ^* is analogous to E _x ^d for the case of delaminations emanating from matrix cracks in the limit the crack spacing D goes to infinity.

$G_{edge}=\frac{1}{2}{\epsilon }_x^{2}h\left(E_{LAM}-E^{*}\right)$

(34)

Figure 8.

Energy release rates for various cases as a function of crack spacing.

As can be seen from Figure 8, cases I and II give higher energy release rate values and would therefore be more likely to occur if there were no defects in the laminate promoting other cases. Case III gives significantly lower release rate values. Case IV is close to Cases I and II and could be favoured if there were inconsistencies in the laminate such as local waviness, porosity etc. It is also seen from Figure 8 that, for relatively small crack spacings, lower than 0.5 mm in Figure 8, Cases I through IV have G values higher than G for edge delamination in equation (34).

This means that for low values of D, the tendency would be for delaminations to start at the tips of matrix cracks. For higher D values, where there are relatively long portions of the laminate without any matrix cracks, the tendency would be to start edge delaminations in those regions.

At this point, the strain energy release rates for a laminate with matrix cracks only and a laminate with matrix cracks and delaminations can be combined to predict delamination onset. This is shown in Figure 9 for a cross-ply laminate. Starting from a laminate under load with matrix cracks, if the load is further increased, energy can be dissipated by forming a new crack or forming a delamination. The point of intersection of the two energy release rate curves in Figure 9 defines two regions of interest: To the left of the point of intersection, the energy release rate for creating a delamination is higher. This means for crack spacings lower than the D _onset value corresponding to the intersection, (D _onset ∼ 0.8 mm in Figure 9), delaminations are the preferred way of energy dissipation. To the right of D _onset , i.e. for larger crack spacings, it is more energy-efficient for the laminate to dissipate energy by increasing the area of matrix cracks. Equating the two strain energy release rates G ^c and G ^c+d gives the value of D _onset . It can be shown that, for a cross-ply laminate, the applied strain which would cause delaminations to start is given by:

${\epsilon }_{a,onset}=\frac{{\sigma }_{x,onset}}{E_x}=\frac{Y_{is}^t}{E_2}\sqrt{\frac{0.27024}{f\left(D_{onset}\right)}}$

(35)where f(D) is given in the Appendix. For a general symmetric laminate under tension and/or shear, the process is basically the same as that just described for cross-ply laminates. However, the solution now cannot be obtained without iterations due to the fact that the equations cannot be normalized with the strains applied to the laminate. This is mainly because of the dependence in the shear modulus of the local applied shear strain, as described in equation (3). As before, the crack spacing D when the two strain energy release rate curves in Figure 9 intersect, is determined through solving:

$G^c\left(N_x,\frac{N_x}{N_{xy}},D\right)=G_{onset}^{c+d}\left(N_x,\frac{N_x}{N_{xy}},D\right)$

(36)where N _x /N _xy is the ratio of applied axial and shear forces on the laminate. But the value obtained from solving equation (36) must also correspond to the matrix crack spacing that the applied load would generate in the laminate. Therefore, equation (7) must also be satisfied. If this is not the case, the applied load combination must be changed until the value of D = D _onset which satisfies them simultaneously is found. Owing to the fact that there are three independent possibilities for G ^c+d _onset from equations (33a-d) out of the four possible delamination configurations examined, there will be a range for D _onset as is shown in Figure 10. Accordingly, three different delamination onset loads can be determined corresponding to the three distinct D _onset values just found. This will give a range of load predictions with the smallest D _onset corresponding to the highest delamination onset load. This would be case III in Figure 10.

Figure 9.

Determination of matrix crack spacing to cause delamination onset.

Figure 10.

Determination of crack spacing at delamination onset.

As there are no openly available data for general laminates, the accuracy of the resulting predictions is assessed by using data on cross-ply laminates. This is shown in Figures 11–13. It should be emphasized here that the predictions have no need of a critical energy release rate value as they are based on critical energy density determined from equation (5) through (7) and (33a-d) making use of Y ^t _is . In Figures 11 and 13, the value of Y^t_is was obtained from Ref. 19:

$Y_{is}^t=1.12\sqrt{2}{Y}^t$

(37)where Y^t is the transverse strength of the material with the Y^t value of 71 MPa in Figure 11 and 58.6 MPa for Figure 13 taken from the references. For Figure 12, the Y^t_is value of 60 MPa was used as reported in the respective reference.

Figure 11.

Onset of delamination load in cracked cross-ply laminates compared to tests from Refs. ² .

Figure 12.

Onset of delamination load in cracked [±25/90_n]s laminates compared to tests from Refs. ²⁰ .

Figure 13.

Predicted crack density at delamination onset for various laminates compared to predictions in Refs. ¹⁴ .

The delamination onset stress for [0/90₂]s and [0/90₄]s laminates is shown in Figure 11. The predictions from the present method are shown next to the test results with a range corresponding to the three Cases examined also shown. The predictions overlap most of the tests showing good agreement. The differences from the test results in Figure 11 are attributed to (a) the equation (37) used for Y^t_is for both laminates as insufficient data were available to capture the dependence of the transverse in situ strength for the lower stack thickness of the 90^o plies, and (b) the delamination pattern reported in Ref. 2 was a mixture of the individual Cases examined here.

Similar good agreement in capturing the trend is observed for [±25/90_n]s laminates when compared to test results from Ref. 20. The comparison is shown in Figure 12. No data scatter were reported in Ref. 20 so no direct comparison of the ranges of the tests to the predictions can be made. As in the previous Figure, Case III corresponds to the highest predicted delamination onset load and Case I to the lowest.

A final comparison is made in Figure 13. The onset of delamination crack density D _onset as predicted by our present method is compared to the predictions by Nairn and Hu ¹⁴ for different values of n in [0_m/90_n]s and [±45/90_n]s laminates as the values of n and m change. Good agreement between the methods is observed.

The efficiency of the presented method allows trade studies to examine the effect of different parameters. One example is shown in Figure 14 where the layup of a laminate under tension is varied. For [0/θ₂]s laminates, values of θ closer to 0^o tend to have higher energy release rate values at delamination onset than values closer to 90^o. This trend however is not monotonic. There is significant stretching/shearing coupling in some of the laminates and the interaction between transverse tension and shear stresses in the cracked θ plies is too complex to lead to a continuously monotonic behaviour.

Figure 14.

Dependence of strain energy release rate at onset of delamination on matrix crack spacing D for different layups.

Other trends that were observed were that the onset of delamination energy release rate in a cracked ply increases with the thickness of the cracked ply or with the thickness of the surrounding uncracked ply and decreases with an increase in the ratio E₁/E₂ of the material.

Delamination growth

If the applied load is increased beyond the point at which delaminations appeared at the tips of matrix cracks, these delaminations will grow. An approach to determine the delamination size as a function of applied load is proposed here. First, it is assumed that the matrix crack spacing does not change when the load increases further. Essentially, it is more efficient for the laminate to dissipate energy in delamination formation rather than creation of additional matrix cracks. This corresponds to crack saturation where, as the loads increases, the spacing remains constant and uniform.^14,23–25

To proceed, a relationship between strain energy release rate for growth and the delamination size is needed. As such a resistance curve, a relatively simple semi-empirical expression proposed by Orange ²⁶ is used:

$G_R=\frac{S^2\pi a}{g\left(1+C_m^{2}a\right)}$

(38)where G _R is the strain energy change per unit area of crack surface created, S and g are the shear strength and modulus of the material, a is the delamination size and C _m is a material-dependent parameter with units 1/√(length). It can be seen that this expression has the standard shape of resistance curve rising sharply and then asymptotically reaching a value for large values of a. It is proposed here to determine the value of C _m by requiring that, for large values of Δa, G _R equals the mode II critical energy release rate G _IIC . Then,

$C_m=\sqrt{\frac{S^2\pi }{{gG}_{IIc}}}$

(39)

The delamination size a _o at delamination onset is determined by equating the strain energy release rate for a laminate with cracks and delaminations, equations (33a-d) to that for a laminate with delaminations only, equation (38):

$G^{c+d}\left(N_{x,onset},D_{onset},a_0\right)=G_R\left(a_0\right)$

(40)

Substituting:

$\frac{1}{4}h\frac{L}{D_{onset}}\left\{\frac{E_x^c\left(\frac{E_x^c}{E_x^d}-1\right)}{{\left[1-\frac{a_0}{D_{onset}}+\frac{a_0}{D_{onset}}\frac{E_x^c}{E_x^d}\right]}^2}{\epsilon }_{x,onset}^2\right\}=\frac{S^2\pi a_0}{g\left(1+C_m^2{a}_0\right)}$

(41)

An accurate approximation for a _o can be obtained by expanding the right hand side in a Taylor series and keeping terms up to second order in a _o . This leads to:

$a_0=\frac{A{C}_m^2-2\frac{A}{D}\left(∆E\right)-1+\sqrt{{\left(\frac{A}{D}\right)}^2\left(2\left(∆E\right)+C_m^{2}D\right)-2\frac{A}{D}\left(C_m^{2}D-2\left(∆E\right)\right)}}{4\frac{A}{D}\left(∆E\right)C_m^2}$

(42)

with C_m given by equation (39) and

$A=\frac{1}{4}\frac{h\hspace{0.17em}g}{S^2\pi }{E}_x^c\left(∆E\right)\left(\frac{L}{D}\right){\epsilon }_{x,onset}^2$

(43)

$∆E=\frac{E_x^c}{E_x^d}-1$

(44)

For delamination growth, equation (40) is solved again but now N _x > N _x,onset replaces N _x,onset and a replaces a _o . As before, there are four cases corresponding to the four cases shown in Figure 4:

$\begin{array}{l}I:a\left(D_{onset},{\epsilon }_x\right)=\frac{A{C}_m^2-2\frac{A}{D_{onset}}∆E-1+\sqrt{{\left(\frac{A}{D_{onset}}\right)}^2\left(2∆E+C_m^2{D}_{onset}\right)-2\frac{A}{D_{onset}}\left(C_m^2{D}_{onset}-2∆E\right)}}{4\frac{A}{D_{onset}}∆E{C}_m^2}\\ II:a\left(D_{onset},{\epsilon }_x\right)=\frac{A{C}_m^2-2\frac{A}{D_{onset}}∆E-1+\sqrt{{\left(\frac{A}{D_{onset}}\right)}^2\left(2∆E+C_m^2{D}_{onset}\right)-2\frac{A}{D_{onset}}\left(C_m^2{D}_{onset}-2∆E\right)}}{4\frac{A}{D_{onset}}∆E{C}_m^2}\\ III:a\left(D_{onset},{\epsilon }_x\right)=\frac{2A{C}_m^2-\frac{A}{D_{onset}}∆E-2+\sqrt{{\left(\frac{A}{D_{onset}}\right)}^2\left(∆E+{2C}_m^2{D}_{onset}\right)-4\frac{A}{D_{onset}}\left(2C_m^2{D}_{onset}-∆E\right)}}{2\frac{A}{D_{onset}}∆E{C}_m^2}\\ IV:a\left(D_{onset},{\epsilon }_x\right)=\frac{2A{C}_m^2-4\frac{A}{D_{onset}}{∆E}^{*}-1+\sqrt{{\left(\frac{2A}{D_{onset}}\right)}^2\left(2{∆E}^{*}+C_m^2{D}_{onset}\right)-4\frac{A}{D_{onset}}\left(C_m^2{D}_{onset}-2{∆E}^{*}\right)}}{8\frac{A}{D_{onset}}{∆E}^{*}{C}_m^2}\end{array}$

(45a-d)

with:

$A=\frac{1}{4}\frac{h\hspace{0.17em}G}{S^2\pi }{E}_x^c∆E\left(\frac{L}{D}\right){\epsilon }_x^2$

(46)

$∆E=\frac{E_x^c}{E_x^d}-1$

(47)

${∆E}^{*}=\frac{E_x^c}{E_x^{del*}}-1$

(48)

Published data with delamination size as a function of applied load are relatively limited. Here, the accuracy of equations (45a-d) is evaluated against test results found in Ref. 3. For the material in Ref. 3, the following values were used: S = 50 MPa, g = 2.85 GPa, and G_IIc = 820 J/m². The comparison of predictions to test results for three cross-ply laminates is shown in Figure 15. As before, a range of predictions is given in Figure 15, corresponding to the four cases of Figure 4. For all three laminates, the tests show there is a region with measured delamination sizes without predictions from the present method. This is because the predictions do not include edge delaminations which would start irrespective of the presence and spacing of matrix cracks. Beyond a certain applied strain, the analysis presented here predicts delamination sizes which are in good agreement with the measured values for the [0/90₄]s and [0/90₆]s laminates and slightly shifted to the left for the [0/90₂]s laminate suggesting that the transverse in situ strength value used in the analysis needs refinement.

Figure 15.

Size of delaminations emanating from matrix cracks compared to tests from Refs. ³ .

Complete history of matrix cracks and delaminations

One of the advantages of the present method is that, being essentially in closed form, it allows very rapid construction of the complete history of cracks and delaminations in a laminate. An example is shown in Figure 16 for the three laminates of Figure 15 under uniaxial tension. In this Figure, crack spacing is measured on the left y axis and delamination length on the right y axis.

Figure 16.

Matrix crack evolution and delamination growth in cross-ply laminates.

For low applied stress, there are no matrix cracks. This is the region with the horizontal lines in the Figure. As the applied stress increases cracks appear and their spacing rapidly decreases. This is shown by the (near) vertical lines in the Figure. Upon further increase of the applied stress, the crack spacing has become uniform and increases more slowly as shown by the curved portion of the three continuous curves in Figure 16. Then, a stress is reached where the cracks saturate and delaminations appear. Their size at onset is smaller than the current crack spacing. The evolution of these delaminations as the stress increases further is delineated by the three dashed curves with the delamination sizes read off the right y axis. At some point, the delaminations grow to the next matrix crack and there is no growth of delaminations beyond that point. This is indicated by the horizontal dashed portions of the three lines.

Discussion

There were several simplifications introduced in the approach presented above which should not be overlooked and which are the subject of future expansions of the work: (a) Thermal (and swelling) stresses were neglected. This was done as a simpler first step to examine the accuracy of the approach without them. Conceptually, the extension to include these stresses is relatively simple by adding the corresponding thermal and swelling strains. It should be noted that, while such an extension is necessary, previous work, ²⁷ has shown that at least for certain laminates and materials, their effect is small. (b) There is a three dimensional stress field at the tips of the matrix cracks and around delaminations whose effects were neglected. Moreover, the fact that, under load the delaminations open the tips of the matrix cracks affecting the local stiffness significantly, was also neglected. These effects can be taken into account with detailed finite element models and are addressed in detail in the work by Kahla. ²⁸ The approach presented here aims at demonstrating that the crack density and the delamination length can be used to predict laminate stiffness in presence of damage or, conversely, that laminate stiffness can be used to determine energy release rate and, through that, delamination onset and growth. The favourable comparison with tests suggests that the method is promising but further improvements are necessary. (c) The three-dimensional stress fields created near the crack tips and delaminations will cause local bending which would be more pronounced if the delamination pattern is unsymmetric as in Case IV. These are neglected in the present work.

Summary and conclusions

The strain energy release rate for symmetric laminates under tension or shear containing matrix cracks or matrix cracks and delaminations was determined. Four different cases were distinguished for the delamination pattern leading to a range of predictions for the crack spacing at which delaminations will start and the corresponding size of delamination onset. A delamination resistance curve was proposed which was used to determine the delamination size for higher applied loads. Good agreement with published test results was observed. It was found that, as the thickness of the cracked ply (or plies) increases, the tendency for matrix crack creation increases and delaminations at the tips of these cracks will start earlier. Also, as the orthotropy ratio of the pristine material increases, the tendency to crack the matrix increases and the load at which delaminations will appear increases. Finally, for a ply under combined shear and transverse tension, the relative magnitudes of these stresses will dictate when cracks and delaminations will appear. The approach presented allows determination of the complete history of matrix cracks and crack tip delaminations for symmetric laminates.