High-fidelity simulations of multiple fracture processes in a laminated composite in tension

Journal of the Mechanics and Physics of Solids 59 (2011) 1355–1373

Contents lists available at ScienceDirect

Journal of the Mechanics and Physics of Solids journal homepage: www.elsevier.com/locate/jmps

High-fidelity simulations of multiple fracture processes in a laminated composite in tension X.J. Fang a, Z.Q. Zhou a, B.N. Cox b, Q.D. Yang a,n a b

Department of Mechanical and Aerospace Engineering, University of Miami, Coral Gables, FL 33124, USA Teledyne Scientific Co LLC, Thousand Oaks, CA 91360, USA

a r t i c l e i n f o

abstract

Article history: Received 14 June 2010 Received in revised form 8 March 2011 Accepted 3 April 2011 Available online 14 April 2011

The augmented finite element method (A-FEM) is used to study the fundamental composite failure problem of delamination and associated damage events spreading from a stress concentrator during tensile loading. The solution exploits the ability of A-FEM to account for coupled multiple crack types that are not predetermined in shape or number. The nonlinear processes of each fracture mode are represented by a cohesive model, which provides a unified description of crack initiation and propagation and can also describe crack coalescence and bifurcation. The study problem is an orthogonal double-notched tension specimen, in which delaminations interact with transverse ply cracks, intra-ply splitting cracks, non-localized fine-scale matrix shear deformation, and fiber breaks. Cohesive laws and constitutive laws for matrix shear deformation are calibrated using literature data from independent tests. The calibrated simulations are mesh independent and correctly reproduce all qualitative aspects of the coupled damage evolution processes. They also correctly predict delamination sizes and shapes, the density of transverse ply cracks, the growth rate of splitting cracks, softening of the global stress–strain curve, and the ultimate strength. A sensitivity analysis relates variability in cohesive law parameters to predicted deviance in engineering properties. Given the known variability in cohesive law parameters, the predicted deviance in ultimate strength agrees with that in experimental data. The importance of including the interactions between different crack systems and nonlocalized shear deformation is demonstrated by suppressing the presence of separate mechanisms; the predicted delamination shapes, splitting crack growth rate, and the stress–displacement relationship fall into significant error. & 2011 Elsevier Ltd. All rights reserved.

Keywords: Fracture Fiber-reinforced composite material Finite elements

1. Introduction Despite extensive research on composite materials, predicting their progressive failure remains a challenging task due to the complexity of the interactions among multiple damage processes. For laminated composites, the damage processes include transverse matrix cracking, fiber rupture (in tension) or kinking (in compression), splitting between fiber and matrix; interlaminar delamination and fine-scale nonlinear shear deformation (Yang and Cox, 2005). Experimentally it has been widely observed that the intra-ply and inter-ply damage modes interact with each other to form complex 3D crack networks. Fig. 1(a) shows the multiple damage modes in a double-notched tension specimen with a symmetric [0/90]s ply

n

Corresponding author. Tel.: þ1 305 284 3221; fax: þ1 305 284 2580. E-mail address: [email protected] (Q.D. Yang).

0022-5096/$ - see front matter & 2011 Elsevier Ltd. All rights reserved. doi:10.1016/j.jmps.2011.04.007

1356

X.J. Fang et al. / J. Mech. Phys. Solids 59 (2011) 1355–1373

Front of delamination

major splitting cracks (H-cracks) in 0°-ply delamination cracks

Splitting cracks in 00-ply Multiple transverse cracking

load axis

spl itt In± ing cr 45° ack plie s s load axis

Fig. 1. (a, b) X-ray radiography reveals damage mechanisms viewed through the ply stack in (a) a double-notch tension specimen with symmetric orthogonal ply stack [0/90]s (Hallet and Wisnom, 2006b), and (b) a quasi-isotropic laminate [  45/þ 45/90/0]s with a circular open hole under tensile loading (Case and Reifsnider, 1999). (c) polarized light micrograph shows distributed damage in a fibrous stitch that has been sheared by a mode II delamination crack (Cox et al., 2008). (d) Scanning electron micrograph of en echelon cracking between two fibers in a carbon–epoxy composite loaded in off-axis tension (Cox et al., 1994).

stack. Dominant splitting cracks in the 01-ply appear as sharply defined horizontal lines and eventually span the specimen. Transverse cracks proliferate in the 901-ply during load increase. In addition, the major splitting cracks are accompanied by wedge-shaped delaminations between the plies (areas of shadow around the splitting cracks). Fig. 1(b) shows the multiple cracking features in a quasi-isotropic laminate ([ 45/þ 45/90/0]s) with a central open hole. The splitting cracks are shorter and the delaminations are lobe-shaped and transverse cracking occurs predominantly in the 7451 plies. The intra-ply cracks, i.e., the splitting cracks in the 901-ply and off-axis cracks in 7451 plies are distributed differently in space in different plies, and show evidence of coupling with delaminations at inter-ply interfaces (in the dark-lobed regions). As well as the crack systems visible in Fig. 1, fine-scale shear deformation has also been observed occasionally within individual plies by microscopic examination (Spearing, 2010). This deformation arises under shear loading. It is continuously distributed, rather than discrete, over spatial scales much smaller than the spacing of any of the crack systems shown in Fig. 1 (the most closely spaced being the transverse ply cracks, with spacing 1 mm to order of magnitude) and is therefore a distinct mechanism. While not always examined in full detail, the distributed shear deformation is likely to consist at the fiber scale of fibrillar crazing or arrays of microcracks, whose lengths typically span the distance between neighboring fibers, i.e., 1 mm to order of magnitude. Examples appear in old experimental studies: Fig. 1c shows fine-scale damage distributed throughout the section of a fibrous reinforcing stitch that has been sheared by a mode II delamination crack; Fig. 1d shows an array of microcracks between two off-axis fibers in a carbon/epoxy composite loaded in tension. In both cases, shear damage arises from mm-scale phenomena. Several authors have concluded that including shear nonlinearity of the ply material is necessary to achieve good correlation between simulations and experiments, especially for predicting splitting crack growth accurately (Wisnom and Chang, 2000; Cox and Yang, 2006; Hallett and Wisnom, 2006a). Thus the total damage system in a polymer composite subject to arbitrary stress states comprises multiple crack types and possibly fine-scale distributed shear damage. Accounting for the strong interactions among these multiple damage mechanisms in an accurate and computationally efficient scheme remains a difficult task. Traditionally, the intra-ply and inter-ply damage processes have been treated separately with different theories. Delamination problems have been analyzed using either linear elastic fracture mechanics methods (e.g., (Pagano and Schoeppner, 2000; Galessgen et al., 2002; Tay, 2003)) or cohesive interface models (Shahwan and Waas, 1997; Wisnom

X.J. Fang et al. / J. Mech. Phys. Solids 59 (2011) 1355–1373

1357

and Chang, 2000; Yang and Cox, 2005; Turon et al., 2006; Xie et al., 2006; Hallett and Wisnom, 2006a). For the intra-ply damage modes, numerous strength based criteria coupled with continuum damage mechanics (CDM) for strength degradation have been used (e.g., (Matzenmiller et al., 1995; Pinho et al., 2006; Lapczyk and Hurtado, 2007; Laurin et al., 2007; Maimi et al., 2007)). Direct coupling of CDM with discrete fracture mechanics models for delamination has also been attempted by several research groups but with only limited success (Choi and Chang, 1992; Carpinteri and Ferro, 2003; Cox and Yang, 2006; Van de Meer and Sluys, 2009a). Moreover, several recent studies have shown that the homogenization process in continuum damage mechanics leads to the loss of key information on multiple-damage coupling at the macroscopic scale and may result in inaccurate prediction of the crack path (Talreja, 2006; Van de Meer and Sluys, 2009a). Unexpected and severe ‘‘stress locking’’ (significant stress remaining in an element that has been assigned complete failure) may also occur (Iarve et al., 2005). The inadequacy of continuum damage modeling has led to the recent trend to integrate explicit representations of all major cracking events into global composite structure models to achieve direct coupling (Iarve et al., 2005; Gonzalez and LLorca, 2006; Hallett and Wisnom, 2006a; Yang et al., 2011). Two critical capabilities that enable the direct coupling are (1) improved cohesive zone models (CZMs) for bulk and interface crack problems to achieve unification of crack initiation and propagation, and (2) improved numerical methods that allow for arbitrary crack initiation and propagation in continua. Although CZMs have been widely used in composite fracture analysis, especially in delamination analyses, their implementation has suffered from a major shortcoming: it has required the potential crack path to be known a priori, so that CZM elements could be directly implanted along the path during the original mesh definition. This greatly limits the application of CZM for problems with a priori unknown cracking paths as frequently observed in composite materials (Fig. 1).

1.1. Methods that allow arbitrary fracture paths Recently several numerical methods have been developed to allow arbitrary cohesive crack initiation and propagation without continuous remeshing. For example, cohesive models have been integrated into the extended finite element method (X-FEM) framework and successfully used to model fracture in homogeneous quasi-brittle materials (e.g., Belytschko et al., 2001; Dolbow et al., 2001; Mo¨es and Belytschko, 2002; Remmers et al., 2003; de Borst et al., 2006). The key feature in the X-FEM formulation is the use of enrichment functions for cracked elements. This is achieved by enhancing the degrees of freedom (DoFs) of all the nodes employed by the elements with internal discontinuity. One of the advantages of this method is that, for an element that encompasses a singular crack-tip, asymptotic displacement functions from fracture mechanics can be used as the enrichment functions and reasonable numerical solutions can be obtained that are partially optimized to the expected asymptotics. One minor problem is that the long-range influence of the singular field is truncated by the use of local element shape functions. As a result, ambiguities arise for those blending elements that connect a tip-enriched element to neighboring elements that are not influenced by the crack. These elements have to be treated specially and sometimes the treatment leads to loss of element locality (Chessa et al., 2003; Fries and Belytschko, 2006). The X-FEM is very effective for treating cracks or even multiple cracks in homogeneous isotropic materials, for which the enrichment functions are known. However, for complex heterogeneous material systems, such as laminated or textile composites, the enrichment functions are not readily available except for some very special cases such as a delamination crack at a symmetric plane in which case the bonded plate/beam can be treated as orthotropic materials and the singular stress field is known (Huynh and Belytschko, 2009). Furthermore, in the cases that multiple cracking systems co-evolve nonlinearly as shown in Fig. 1, X-FEM will have difficulty in dealing with the coupling among the multiple cracks. A promising alternative numerical method for handling arbitrary cracking problems is the phantom node method. This line of development follows the original work of Hansbo and Hansbo (2004), who first established that an arbitrary discontinuity in an element can be accounted for by adding an extra element on top of this element, with each of these elements accounting for the stiffness and force contributions from the bisected physical domains. The two discontinuous domains are connected by linear or nonlinear springs (Hansbo and Hansbo, 2002, 2004), or cohesive failure tractions (Mergheim et al., 2005, 2007). The addition of element is typically realized by adding nodes that are geometrically identical to the original corner nodes (Song et al., 2006; Van de Meer and Sluys, 2009b; Van de Meer et al., 2010). One major advantage of this method is that it uses only standard finite element (FE) shape functions, which renders it completely compatible with standard FE programs. However, adding elements as a crack evolves changes the numerical model dimension dynamically, thus requiring full access to the source code, which is not always possible with commercial codes (the same is also true for X-FEM). To ease this inconvenient restriction, Ling et al. (2009) proposed an element with double nodes that can treat an arbitrary intra-element discontinuity. They mathematically proved that such a physical element, if containing a discontinuity (e.g., a cohesive crack or a bi-material interface) can be augmented internally (hence the name augmented finite element—A-FEM) as two combined mathematical elements (MEs), each using standard FE shape functions only. Further, they demonstrated that material heterogeneity in composite materials can be conveniently considered with this element. The initial success of this method was demonstrated using an open-sourced FE code (Zienkiewicz and Taylor, 2005). In this paper, the A-FEM proposed in Ling et al. (2009) will be adpated into a formulation that can be directly implemented into any commerical codes with limited, or even without, direct access to source code.

1358

X.J. Fang et al. / J. Mech. Phys. Solids 59 (2011) 1355–1373

1.2. Cohesive and bridged crack models of fracture Apart from details affecting the convenience of implementation, the X-FEM and A-FEM differ importantly in their physical representation of fracture: X-FEM assumes the presence of a singular field at the crack-tip while A-FEM assumes that stress fields remain finite and continuous at the crack tip. In the literature, the former closure condition has been referred to as a ‘‘bridged crack’’ model, the latter as a ‘‘cohesive model’’ (Cox and Marshall, 1994). When present, a singular tip field poses special problems for modeling crack initiation and bifurcation. The tip field influences the criterion for crack propagation, but is not present until the crack exists. In a cohesive model, the same cohesive constitutive law governs crack initiation and propagation. In a bridged crack model, initiation is often assumed to have occurred; the material contains some initial defect. If not, then initiation must be determined by some criterion expressed in terms of finite fields, after which a transition must be implemented to a crack of nonzero length with its tip singularity. Executing this transition in a numerical simulation is problematic and can introduce mesh dependence. The cohesive model, by unifying initiation and propagation criteria, greatly simplifies handling the transition from initiation to propagation, crack bifurcation, and the impingement of one crack upon another. These benefits will be demonstrated below. Nevertheless, using a cohesive model of fracture, especially one with a relatively simple cohesive law, and omitting any crack tip singularity cannot be universally valid for all materials. Some materials, such as fiber-reinforced brittle matrix (ceramic) composites, clearly possess distinct tip and wake fracture processes. To get the fracture mechanics right, one must represent both processes (Cox and Marshall, 1994). Either a bridged crack model must be used or a cohesive model with a sharp peak at zero displacement to represent the tip process. These models are physically equivalent in appropriate  1996), but the question of which is more convenient in large-scale computations remains limits (Carpinteri and Massabo, to be resolved. We do not address this important question here, but restrict modeling to relatively simple cohesive laws. Aspects of fracture in polymer composites for which two-part (two-process) cohesive laws might be required are under current study by the authors. 1.3. Outline of paper In Section 2, the element augmentation of Ling et al. (2009) is realized in a formulation that uses a standard finite element description without access to the source code. Implementation into a commercial FEM package (ABAQUS) demonstrates its full compatibility with standard FE packages. In Section 3, mesh independence is established and the crucial role of the cohesive laws in fracture behavior is tested by a sensitivity analysis, varying parameters around the values set by calibration against independent literature data. Section 4 shows that high fidelity is achieved only by including all fracture modes together with non-localized matrix shear deformation. 2. Augmented finite element method 2.1. Augmentation of displacement field for intra-element discontinuity Consider an element that is cut by a cohesive crack as shown in Fig. 2(a). Standard FE cannot describe the discontinuous displacement fields between the two severed domains. In the A-FEM method developed by Ling et al. (2009) with double nodes (one set of physical nodes 1–4, and another set of ghost nodes 10 –40 ) as shown in Fig. 2(a), the two severed physical domains can be separately approximated by two mathematic elements (MEs) as shown in Fig. 2(b) and (c). The active material domain in each ME is indicated by the shaded area. Stiffness and force integration is performed only over the associated active domain in each ME. The two MEs are then connected by a cohesive failure

4 (4’)

3 (3’) Ω2

Γc2

1 (1’)

t1

(7) (6)

Γc1e

Discontinuity

3

4 Γ2e

e

(8) (5)

3’

4’

e

Ω1

e Γ c1

Ω2

e

Ω2

e

Γ c2

-t1

Ω1 e

e

Γ1 e 2 (2’)

1

ME 1

2

1’

2’ ME 2

Fig. 2. (a) An A-FEM with double nodes traversed by an intra-element cohesive crack. This element can be treated by defining two mathematical elements (b, c), each having the same geometrical dimension as the original element but with different domains (material volumes) for stiffness integration (Ling et al., 2009).

X.J. Fang et al. / J. Mech. Phys. Solids 59 (2011) 1355–1373

1359

description along the crack line uðxÞ ¼

n

NðxÞf1 ðxÞ

NðxÞf2 ðxÞ

( ) o d1

ð1Þ

d2

where NðxÞ are the standard FE shape functions and da (a ¼1 and 2) are the nodal displacements of ME-a. f1(x) and f2(x) are two functions defined below to ensure only the physical domains in the two MEs are used in stiffness and force integration: ( 1 x 2 Oe1 1þ HðxÞ 1HðxÞ ; f2 ðxÞ ¼ ; HðxÞ ¼ f1 ðxÞ ¼ 1 x 2 Oe2 2 2 Here H(x) is generalized Heaviside step function. Eq. (1) now allows for different displacement fields in physical domain

O1 e and O2 e , which makes it possible to account for the discontinuous displacements (or, displacement jumps) across the cohesive crack Gc. The displacement jumps can be conveniently computed as  wðxÞ ¼ u2 ðxÞu1 ðxÞ ¼ N

N

( )  d1 d2

ðx 2 Gc Þ

ð2Þ

If an element does not host a cohesive crack and is not influenced by a crack, either f1(x) or f2(x) is zero in the associated ME. Then Eq. (1) reduces to the standard FE shape function interpolation. Numerically, this can be simply enforced by equating all the DoFs of the ghost nodes to those of the associated physical nodes. In laminated polymer–matrix composites made of unidirectional plies, an intra-ply crack, being a transverse matrix crack or a splitting crack, always propagates along the local fiber direction in each individual ply. Therefore, it is advantageous to express the displacement jump in the local coordinates as shown in Fig. 3. In this study, the out-of-plane deformation is neglected because the system is dominated by in-plane deformation. Under the local coordinates the displacement components can be obtained using a rotational matrix R, i.e.,  T w0 ¼ dn ds dt ¼ Rw ð3Þ where R is the rotational matrix between the global and local coordinates, i.e., 2 3 sin y cos y 0 6 7 R ¼ 4 cos y sin y 0 5 0 0 1

ð4Þ

Similarly, the cohesive traction along the discontinuity can be expressed in local coordinates as t0 ¼ ftn0 ,ts0 ,tt0 g ¼ t0 ðw0 Þ

ð5Þ

The cohesive law for mixed mode fracture used in this study is given in Appendix A. 2.2. Discretized weak form The weak form of the equilibrium equations for static analysis without body force can be derived from the virtual work principle (Ling et al., 2009), Z ( T ) Z ( NT Ff ðxÞ ) Z ( BT rf ðxÞ ) N t 1 1 dO þ dG dG ¼ ð6Þ T T T Oe B rf2 ðxÞ Gc N t GF N Ff2 ðxÞ where B is the strain matrix, r is the stress matrix, F is the external force array, and tðwÞ ¼ RT t0 ðw0 Þ is the cohesive traction array. The discretized form of Eq. (6) is nonlinear in general because the cohesive traction in the second integral is typically a nonlinear function of the cohesive opening w. However, it is straightforward to linearize the equation using an

t

z

s

y n x

Fig. 3. Local coordinates (n, s, t) defined by the cohesive crack with n the in-plane direction normal to fiber orientation, s the in-plane direction aligned with fiber orientation and t the out-of-plane direction.

1360

X.J. Fang et al. / J. Mech. Phys. Solids 59 (2011) 1355–1373

incremental scheme. The linearized equation can be written as 9 " #( ) 8R < G NT DFf1 ðxÞ dG = K11 K12 Dd1 F ¼ R : G NT DFf2 ðxÞ dG ; K21 K22 Dd2

ð7Þ

F

where Dda (a ¼1, 2) are incremental nodal displacements and DF is the external load increment. The matrices Kij are Z Z 2 K11 ¼ BT RT D0 RBf1 ðxÞ dO þ NT RT D0coh RN dG K22 ¼

Z

Oe

Gc

T

Oe

T

0

B R D

2 RBf2 ðxÞ dO þ

Z

Gc

NT RT D0coh RN dG

Z K12 ¼ K21 ¼  NT RT D0coh RNdG

ð8Þ

Gc

where D0 is the ply stiffness matrix in local coordinates and D0coh is the instantaneous tangential stiffness matrix of the cohesive law, which is given in Appendix A. Eqs. (7) and (8) are fully compatible with any standard FE program including commercial codes with limited access to the source code. In this study, the A-FEM has been implemented in the commercial code ABAQUS (V6.7) as a user-defined element. 3. High-fidelity simulations of a double-notched tension specimen 3.1. Specimen geometry and numerical model In this section, the fidelity of the A-FEM and cohesive fracture models in simulating multiple cracking systems in laminated composites will be assessed through comparison with experimental data due to Hallett and Wisnom (2006a, 2006b). Their specimen geometry and notch configuration are shown in Fig. 4(a). Displacement-controlled in-plane tensile loading was applied. The two dots on the specimen surface acted as reference points for measuring displacements on the specimen itself, to eliminate test machine compliance from the load–displacement data. Due to symmetry in specimen geometry, ply stacking and loading, only a quarter of the specimen is modeled with symmetric boundary conditions applied at the left and bottom edges, as shown in Fig. 4(b). The possibility that damage might develop asymmetrically is therefore not considered; damage in the tests being simulated did remain close to symmetric. Further, in this study since the deformation in the planar DNT specimen under simple tension is mostly within the x–y plane, each ply is modeled by plane stress elements with ply thickness properly assigned. The out-of-plane stresses are therefore neglected. The ply thickness is 0.127 mm (0.005 in.) Transverse ply cracks and splitting cracks were modeled using A-FEM elements within the plies: the flexibility of A-FEM therefore allows these two crack systems to arise at arbitrary locations during simulations. Since the potential fracture plane of the delaminations is known, they were treated differently by pre-locating cohesive zone elements along all ply interfaces. This mixed meshing strategy allows the complex interaction between intra-ply cracks and inter-ply delamination cracks to be represented with appropriate realism. Using pre-planted cohesive elements introduces a risk to fidelity: if they are not correctly located, their use can significantly delay the initiation of cracks, because a pre-planted cohesive element may not be exactly at the location identified for initiation by the local initiation criterion. The A-FEM formulation will always find and use the maximum local stress for initiation (to within the accuracy of the stress representation). This risk was of most concern for the transverse ply and splitting cracks.

w = 20mm

θ = 60° 20.8 mm 2a/w = 0.5

End tabs

a L = 100mm

Fig. 4. (a) Geometry of the double-notched tension specimen tested by Hallett and Wisnom (2006b). The two filled dots in the specimen were used as markers for direct optical measurement of the displacement, so that the machine compliance could be eliminated from experimental data. (b) The quarter geometry numerical model used in our A-FEM simulation.

X.J. Fang et al. / J. Mech. Phys. Solids 59 (2011) 1355–1373

1361

Because the modeling strategy uses plane stress elements for plies, it implicitly depicts any intra-ply crack (transverse or splitting crack) as immediately penetrating the entire ply thickness, once initiated. It therefore starts to interact immediately with interface cohesive elements. This treatment is valid as long as the mode-I toughness (which controls transverse cracking in 901-plies) is smaller than the mode-II toughness (which controls delamination), which is the case in most polymer composites (including the subject composite). A detailed study of a single transverse ply crack and associated delaminations shows that, given this ordering of the toughness values, the intra-ply cracks do span the ply before delamination is triggered (Zhou et al., 2010). This point will be elaborated below.

3.2. Mesh sensitivity Mesh-independent results can only be expected in a nonlinear fracture analysis if the mesh is fine enough to depict the cohesive damage process zone adequately. An adequate mesh refinement can be estimated using the analytical approximations for the cohesive zone size that are summarized by Yang and Cox (2005). For the calibrated cohesive parameters in Table 1 below, the estimated cohesive zone length for mode-I fracture (transverse intra-ply cracks) is lIcoh ¼ 0:16 mm and that for mode-II fracture (e.g., splitting cracks) is lIIcoh ¼ 0:53 mm. Accordingly, mesh dependence was studied using the three mesh resolutions of le ¼0.1, 0.2, and 0.4 mm to represent the numerical domain in Fig. 4(b). Comparisons were made for the predicted stress–extension curve (Fig. 5a), the splitting crack growth vs. applied nominal stress curve (Fig. 5b), and the transverse crack density (Fig. 5c–e). Specimen compliance and splitting crack growth: from Fig. 5(a) and (b), the stress–displacement curves and splitting crack length vs. applied stress curves obtained for all three meshes are very consistent. This is consistent with the previous conclusion that A-FEM results are mesh independent provided that the mesh size is smaller than the cohesive zone size (Ling et al., 2009). Transverse crack spacing: in the present theory, the locations and spacing of transverse intra-ply cracks is determined automatically according to the computed stress state at any time in a simulation; unlike other recent works (e.g., Van de Meer et al. (2010)), no constraint is imposed on crack spacing other than the fact that one element and its nodes can support only one crack. Therefore, another crack cannot propagate into any of the (neighboring) elements that share the nodes of the already cracked element and the limit crack spacing in the mesh of Fig. 5 is that any two transverse cracks must be separated by at least one element. That is, the minimum possible transverse crack spacing is 0.2, 0.4, and 0.8 mm for meshes le ¼0.1, 0.2, and 0.4 mm, respectively. The predicted cracking systems for the three meshes are shown in Fig. 5(c–e). The predicted transverse crack spacing for meshes with le ¼0.1, 0.2, and 0.4 mm are, respectively, 0.9 mm (9le), 1.0 mm (5le), and 1.2 mm (3le). The changes in crack spacing are small given the large variation in mesh resolution and acceptably small in that they are less than variations that might be expected from other factors, such as stochastic local material strength. In addition, the predicted delaminations (dark regions in Fig. 5c–e) are also all very consistent. The demonstrated mesh independence also implies that the numerical formulation presented in Section 2 is able to capture correctly the nonlinear damage coupling among the splitting, transverse, and delamination cracks.

3.3. Nonlinear damage models and their calibration The progressive damage in the DNT specimen was reported in detail by Hallett and Wisnom (2006a, 2006b). There were five main types of damage processes evolving in a coupled manner before final failure occurred (Fig. 4a): (1) splitting cracks in the 01-ply emanating from the notch tip and propagating along the loading direction, which started when the applied nominal stress level reached about 25% of the final strength (260–290 MPa); (2) multiple transverse cracking in the 901-ply with a saturation cracking density of about 2 cracks/mm; (3) triangular delaminations propagating between the 01- and 901-plies in a self-similar fashion with wedge angle of 7–101; (4) shear nonlinearity; and (5) fiber rupture responsible for ultimate failure of the specimen. In addition, data for the load vs. the displacement between the two reference points shown in Fig. 4a (filled dots), and the load vs. the splitting crack length as a function of applied nominal stress, were available for comparison. All these observations and data will be used to test predictions. Each observed mechanism is capable of arising in the simulations, but its presence is not prescribed in advance. Table 1 Transversely isotropic laminar elasticity and baseline cohesive parameters used in this study. Laminar elastic properties

Baseline CZM parameters

E1 (GPa)

E2 (¼E3) (GPa)

n12(¼ n13) (dimensionless)

n23 (dimensionless)

G12(¼ G13) (GPa)

43.9

15.4

0.3

0.3

4.34

Mode-I

Mode-II

Mode III

GIC (J/m2)

s^ n (MPa)

GIIC (J/m2)

t^ s (MPa)

GIIIC (J/m2)

t^ t (MPa)

250

120

900

80

900

80

300 250 200 150 100 50 0

le = 0.2 mm le = 0.1 mm le = 0.4 mm

0

0.1 0.2 0.3 0.4 Displacement, mm

0.5

Length of splitting crack, mm

X.J. Fang et al. / J. Mech. Phys. Solids 59 (2011) 1355–1373

Nominal stress, MPa

1362

12 10 8 6 4

le = 0.4 mm

le = 0.1 mm le = 0.2 mm

2 0 0

300 100 200 Nominal stress, MPa

400

Fig. 5. Comparison of predicted numerical results with three different mesh resolutions: (a) stress–displacement curves, (b) splitting crack length vs. applied nominal stress curves, and (c–e) coupled multiple cracks include transverse cracks (vertical lines), splitting cracks (horizontal lines), and interlaminar delamination (dark triangular zones) at a stress level of 238 MPa.

This section details the constitutive models that were used to represent elasticity and the various damage modes. All parameters in each constitutive model were calibrated using data that were independent of the simulated tests.

3.3.1. Elasticity Each ply is assumed to be homogeneous and transversely isotropic. Calibration: the transversely isotropic ply properties are determined by comparing simulations to experimental measurements of composite stiffness. For the subject materials (E-Glass/913), this calibration was carried out by (Cui and Wisnom, 1992), whose results are shown in Table 1.

3.3.2. Fiber rupture Fiber rupture in 01-plies is introduced by a material degradation model that modifies the affected computational element: ( E1 ðe1 r ef Þ E¼ ð9Þ  0 ðe1 4 ef Þ where E1 is Young0 s modulus along the fiber direction and ef is a threshold strain value along the fiber direction at which fiber rupture occurs. Calibration: the critical strain, ef, is chosen to be 3.5%, a value determined empirically for an E-Glass/913 composite by Hallett and Wisnom (2006b). It is typical of glass fibers.

X.J. Fang et al. / J. Mech. Phys. Solids 59 (2011) 1355–1373

1363

3.3.3. Fine-scale shear deformation A tension test of a 7451 laminate generates shear stress,t12 aligned with the fibers within individual plies (Jelf and Fleck, 1994). The simple relation t12 ¼ sa, where sa is the applied stress, allows a constitutive law for the ply material in shear to be deduced from the global stress–strain data of the test (see also the comparison with Iosipescu shear tests by Kumosa et al. (2002)). Prior research on kink band formation encourages the use of this calibration: the critical shear stress level deduced from 7451 tension data can be used to predict the onset of kink bands in composites in compression. A kink band is a shear-mediated damage mechanism that acts within fiber bundles via the types of damage seen in Fig. 1c and d (Cox et al., 1994; Jelf and Fleck, 1994). The constitutive law calibrated by a 7451 test for the subject material is shown in Fig. 6 (Wisnom, 1994). The observed softening sums the combined effect of fine-scale plasticity of the type shown in Figs. 1c and d as well as discrete delaminations and intraply cracks, whose size and spacing are associated with larger spatial scales. However, the delaminations and intraply cracks contribute to nonlinearity mainly at the higher strain levels in Fig. 6, whereas the simulations presented below (e.g., Fig. 9) predict that plastic strains are small (o0.02) over most of the material. Furthermore, the domain in Fig. 9 that contains larger plastic strains (but stillo0.04) is small in spatial extent compared to the predicted spacing of transverse intraply cracks. Therefore, the mechanism of continuously distributed plasticity represents phenomena that are either occurring at strains beneath those at which large-scale cracks act in the calibration test or are too spatially confined to be included in the effects of transverse cracks. Therefore, including shear plasticity and transverse cracks separately in the simulations does not lead to double counting of a single phenomenon. Consistently, the characteristic of the data of Fig. 6 that dominates the fitting procedure of the calibration is the strain at which the knee occurs (E0.02). The data of Fig. 6 are well approximated by the Ramburg–Osgood law: !n ! t012 jt12 j jt12 j g12 ¼ sgn½t12  0 þ a 0 ð10Þ G12 G12 t012 where G012 is the undamaged in-plane shear modulus and t012 is the limiting elastic shear stress. Calibration: calibration of a constitutive law for distributed nonlinear shear deformation was achieved by fitting Eq. (10) to the data from (Wisnom, 1994). In the fitting function, the undamaged in-plane shear modulus G012 was assigned the value 4.34 GPa and the limiting elastic shear stress t012 (which is related to the yield or flow stress) was assigned the value 30 MPa, using data from (Wisnom 1994). With G012 and t012 chosen thus, fitting Eq. (10) to the data of Fig. 6 yielded a ¼0.01 and n ¼6.8. The stress–strain data are fitted within 2% over the entire test range.

τ12

3.3.4. Cohesive model for delamination, splitting, and transverse intraply cracks The mixed-mode cohesive model detailed in Appendix A was employed to represent the nonlinear process zone in all types of cracks. The law prescribes the traction vector for any crack displacement vector. In this law, the modes are represented by separable functions and identical laws are assumed for modes II and III for the reasons discussed by Yang and Cox (2005). In its most general formulation, the law for each mode is quadrilateral in shape, but in the present work it was reduced to a triangular shape by making the slopes of the two softening zones of Fig. A.1 equal. A criterion that is based on the total energy of deformation specifies the complete failure of the cohesive mechanism (all traction components go simultaneously to zero) for any given mode ratio (Appendix A). In evaluating the cohesive law for delamination cracks, out-of-plane deformation (mode I) is not included. In applying the cohesive law to splitting cracks and transverse ply cracks, which may arise at any location, an initiation criterion must also be stated. The initiation criterion is specified using parameters of the cohesive law, because it refers to the same material physics. Once any crack is initiated, its propagation is also governed by the cohesive law.

90 80 70 60 50 40 30 20 10 0

Experimental data (Wisnom 1994) Ramberg & Osgood fit

0

0.02

0.04

0.06

0.08

0.1

γ12 Fig. 6. Shear stress–strain data from (Wisnom, 1994) fitted by a Ramberg–Osgood law. The shear stress plotted refers to an axis system aligned with the fibers in a 451 ply, and is therefore one-half the applied tensile stress. The strain is the engineering shear strain.

1364

X.J. Fang et al. / J. Mech. Phys. Solids 59 (2011) 1355–1373

In applying the cohesive law to delamination cracks, no initiation criterion need be specified, because the cohesive elements used to model delamination are already present on the known potential fracture plane. The initiation criterion used for splitting and transverse cracks accounts for the effects of tension normal to the fracture plane and in-plane shear, viz.,    2 hsn i 2 ts þ ¼1 ð11Þ s^ n t^ s where sn and ts are the normal tensile stress and the in-plane shear stress, averaged over an element, with the coordinates n perpendicular to local fiber direction, s along the fiber direction, and both in the plane of the laminate. The Macauley bracket hi ¼ max½,0 appears so that normal compression (sn o 0) cannot induce crack initiation. The form of Eq. (11) reflects the fact that, in laminated polymer composites, splitting and transverse cracks tend to follow the local fiber ^ n and t^ s , are the direction (Fig. 1) and propagate under constrained, mixed-mode conditions. The material parameters, s maximum values of modes I and II tractions, respectively, in the calibrated cohesive law. Calibration: cohesive laws and thus the parameters appearing in the initiation criterion were calibrated by using published data to specify the material toughness (the area under the cohesive law) for each of mode I and mode II and the ^ n and t^ s (Eq. (11) and Fig. A.1). The mode III law was set to be identical to the mode II law. Given the material strengths s ^ n or t^ s , the maximum displacement for nonzero tractions is also area under the curve and the maximum traction s determined. The initial steep linear increasing portion of the law (Fig. A.1) is included to control numerical instabilities and has no significance for the fracture process. Identical calibration parameters were assumed for all three crack types: delamination cracks, splitting cracks, and transverse ply cracks. The mode-I toughness of the ply material (E-Glass/913) was determined empirically using tensile fracture tests on single plies with the load axis perpendicular to the fiber direction, yielding GIC ¼250 J/m2 (Hallett and Wisnom, 2006a). The mode II toughness for a laminate made from the same materials was determined by analyzing end-notch flexure delamination tests, yielding two reported values GIIC ¼900 and 1040 J/m2 (Hallett and Wisnom, 2006a, 2006b). The baseline calibration values were taken to be GIC ¼250 and GIIC ¼900 J/m2, with expected variance of 10%. ^ n , primarily influences transverse ply cracks, since both splitting and In the subject problem, the mode I strength, s delamination cracks are shear dominated. In the lay-up studied, the transverse 901-plies are surface plies constrained only on the side that is bonded to the interior 01-plies. Therefore, transverse ply cracks appear as surface channeling cracks, whose behavior has been studied by both the composites and thin film communities (Dvorak and Laws, 1987; Thouless, 1990; Davila et al., 2005; Camanho et al., 2006). According to Camanho et al. (2006), the local stress in the transverse ply required to initiate a surface channeling crack is approximately sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 4GIC ^ ð12Þ sn ¼ phð1=E2 n221 =E1 Þ ^ n from the calibrated where h is the ply thickness (0.127 mm). Since the elastic constants are known, Eq. (12) predicts s ^ n ¼141 MPa. However, this value should be regarded as an upper-bound estimate, value of the mode I toughness, yielding s because Eq. (12) assumes a perfectly bonded interface between neighboring plies. Furthermore, the transverse strength of unidirectional plies of E-Glass/913 exhibits specimen-to-specimen discrepancy of  15% (Wisnom and Jones, 1996). Summing these considerations, the value 120 MPa was used as a baseline value in the initiation criterion, Eq. (11), and in the mode I cohesive law, with expected variance of 15%. This same calibration was also used for the mode I component of the cohesive zone in splitting cracks. Since the splitting crack is an embedded crack (in the central 01 plies), its configuration is a tunneling crack, rather than the channeling crack to which Eq. (12) refers. However, in the particular lay-up studied here, the 01 plies are exactly twice the thickness of the 901 plies; and the effect of the different thickness on the expected critical stress exactly cancels the difference expected due to the fact that the splitting crack tunnels rather than channels (Camanho et al., 2006). Therefore, use of the same mode I calibration is consistent. Prior estimates of the shear cohesive strength, t^ s , for the subject E-Glass/913 material have been based on the shear deformation data of Fig. 5, with the value 75 MPa picked out by Hallett and Wisnom (2006a)) and 70 MPa by Wisnom and Chang (2000). The uncertainty in the appropriate value reflects the range of the shear stress beyond the value at the knee of the data (Fig. 6). In this work, t^ n was assigned the baseline calibration value of 80 MPa, with expected uncertainty of 10%. The same value was used for all crack types. The set of baseline calibration values for cohesive law parameters is shown in Table 1. Variation of the parameters from these values will be made to assess sensitivities. 3.4. Numerical prediction and validation The experimental images taken at three different load levels show the evolution of three major cracking systems (Fig. 7(a)). The damage evolution predicted at the same load levels using the baseline calibration values of Table 1 for the cohesive parameters is shown in Fig. 7(b). In this figure, the splitting crack is shown as a solid white line, the transverse cracks as fine white or red lines, and the delamination by the black zone attached to the splitting crack. The A-FEM simulation successfully

X.J. Fang et al. / J. Mech. Phys. Solids 59 (2011) 1355–1373

1365

Nominal Stress, MPa

400

No transverse cracking Single transverse cracking

Typical exp. curve

300 Arbitrary transverse cracking

200

Experimental strength range (Hallet & Wisnom 2006)

εf = 3.5%

100

ˆ = 120 MPa; GIc = 250 J/m2 ˆ = 80 MPa; GIIc = 900 J/m2

0 0

0.2 0.4 Displacement, mm

0.6

Splitting crack length (mm)

Fig. 7. (a) Three snap shots of the DNT specimen under different load levels (25%, 65%, and 100%), showing the evolution of the coupled transverse cracks (vertical lines), splitting cracks (horizontal line) and delamination zones (triangular dark zones). (b) the A-FEM predicted crack evolution at approximately the same load levels.

12

No transverse cracking εf = 3.5% ˆ = 120 MPa; GIc = 250 J/m2 ˆ = 80 MPa; GIIc = 900 J/m2

8

Single transverse cracking

Experiment (Hallet & Wisnom 2006)

4

Arbitrary transverse cracking Hallet & Wisnom (2006)

0 0

100 200 Nominal Stress (MPa)

300

Fig. 8. Comparison of experimental results reported by Hallett and Wisnom (2006b) with predictions: (a) nominal stress vs. displacement curves for three cases: (1) with arbitrary transverse cracks, (2) with single transverse cracks at notch tip; and (3) no transverse crack. (b) Splitting crack length as functions of nominal stress for the same three cases. The numerical results obtained by Hallett and Wisnom (2006a) using pre-planted CZM elements at the notch tip are also shown. Note: The experimental curve in Fig. 8(a) and the following Figs. 10–13(a) is a re-claibrated one using the two reference dots in Fig. 7(a), kindly provided by Dr. S.R. Hallet.

reproduces all the important features of the three cracking systems. The self-similar delamination zone angle is predicted to be 81, in agreement with the experimental value of 7–101. The splitting and multiple transverse cracks are also well captured. The predicted transverse crack density is  1 crack/mm, which is about half of the measured density. This can be explained as a consequence of assuming symmetry in the damage in simulations. In experiments, transverse cracks are spawned from each of two splitting cracks, so that the final crack density in the region between the two splitting cracks can be up to twice the initiation density along either crack. This possible doubling effect is absent from the simulations, because only one quarter of the specimen, containing only one splitting crack, was simulated. 3.4.1. Load and splitting crack length vs. applied displacement The simulations predict the entire stress–displacement curve up to final failure, and the splitting crack length as a function of applied nominal stress, to within anticipated material variance. The predicted stress as a function of the displacement is compared directly with experimental data in Fig. 8(a). The simulation with arbitrary cracking reproduces the experimental curve very well: it follows the entire experimental curve, including the mild nonlinearity beyond the initial elastic regime, and, for this particular specimen, predicts the ultimate strength quite accurately. The splitting crack growth as a function of applied stress is shown in Fig. 8(b). In this figure, each of the experimental data points is an average of five tests reported by Hallett and Wisnom (2006b). The predicted splitting crack growth curve with all mechanisms enabled in the simulations again follows closely the experimental data during the entire loading stage. 3.5. Investigation of damage coupling, sensitivity to cohesive parameters, and influence of transverse ply cracks To show the necessity of including transverse cracks in strength analyses, predictions of stress vs. strain are also shown in Fig. 8a for two modified cases, viz., (1) without considering transverse cracking, and (2) with only one transverse crack

X.J. Fang et al. / J. Mech. Phys. Solids 59 (2011) 1355–1373

notch

1366

12

x2

0.04 x1

0.02 0.0

Nominal Stress, MPa

400 experiment (Hallet & Wisnom 2006b) Linear shear

300

nonlinear shear

200 100

f = 3.5% ˆ = 120 MPa; GIc = 250 J/m2 ˆ = 80 MPa; GIIc = 900 J/m2

0 0

0.2 0.4 Displacement, mm

0.6

Splitting crack length (mm)

Fig. 9. Snapshot of distributed fine-scale shear deformation in the 901 ply. The strains are meaningful ahead of the tip of the splitting crack (horizontal white line). Two transverse cracks (vertical white lines—the one at notch tip is not visible) are present at this load level. Similar distributions are also predicted in the 01-ply. The large apparent shear strains (gray zones) immediately adjacent to the splitting crack are an artifact of a plotting deficiency of ABAQUS and should be ignored.

12

Experiment (Hallet f = 3.5% & Wisnom 2006) ˆ = 120 MPa; GIc = 250 J/m2 ˆ = 80 MPa; GIIc = 900 J/m2 Nonlinear shea r

8

Linear shear

4

0 0

100 200 Nominal Stress (MPa)

300

Fig. 10. Influence of shear nonlinearity on (a) the stress vs. displacement relation and (b) splitting crack growth. (c) the A-FEM predicted crack evolution without shear nonlinearity at three load levels. The splitting crack shown as white solid line. The transverse cracks are indicated by the white or red colored lines; while the delamination is shown by the black zone attached to the splitting crack (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article).

permitted (a simplification considered by Wisnom and Chang (2000) and Hallett and Wisnom (2006a)). Whereas the prediction with arbitrary transverse cracking is accurate, if no or only one transverse crack is allowed the global stiffness reduction beyond the proportional limit cannot be adequately captured. In contrast, the propagation of the splitting crack is not very sensitive to the transverse crack details, with similar predictions for arbitrary transverse cracking, a single transverse crack, or no transverse cracks (Fig. 8b). This perhaps explains why the simulation by Hallett and Wisnom (2006a), in which a single cohesive crack was pre-planted at the notch tip, also did a reasonable job of capturing the overall splitting crack growth. However, if cohesive elements are not placed accurately at the location of maximum stress, splitting crack initiation can be significantly delayed (Fig. 8b). Further simulations in the present study confirmed the sensitivity of splitting crack initiation to assumed location. (No location was assumed in the preferred, high-fidelity formulation and simulations.) 3.5.1. Influence of shear nonlinearity Fine-scale distributed shear deformation influences all fracture modes. The shear strain distributions in the surface 901-ply are shown in Fig. 9, with splitting and transverse cracks superimposed. The largest shear strains are concentrated in a narrow region that extends ahead of the splitting crack tip. As demonstrated by the single transverse crack seen in Fig. 9, which has just initiated and is incompletely propagated, any transverse crack is initiated in material that is already damaged by distributed shear deformation. Therefore the transverse cracks are initiated under mixed-mode conditions, rather than the approximately mode-I conditions that would exist if they initiated in undamaged material far from the splitting crack.

X.J. Fang et al. / J. Mech. Phys. Solids 59 (2011) 1355–1373

1367

When distributed shear deformation is omitted from simulations, the predicted damage pattern is qualitatively wrong. Fig. 10(c) shows the coupled damage evolution for three load levels similar to those in Fig. 7(b) without shear nonlinearity. The delamination zone is much larger, fewer transverse cracks have developed, and they lag significantly behind the splitting crack tip. Without shear nonlinearity, energy release is directed more towards delamination, rather than to the initiation of transverse cracks. The extensive delamination results in less energy being available for driving splitting crack growth, and thus slower splitting crack growth at large stress levels. Similar numerical observations were reported in (Wisnom and Chang, 2000) for another material system (T300/914), although in that study transverse cracking was not considered. Shear nonlinearity has a correspondingly strong influence on the stress–displacement curve (Fig. 10). Comparing the experimental data with the stress–displacement curves obtained with the calibrated shear nonlinearity included or excluded, one may conclude that the initial specimen stiffness is not very sensitive to shear nonlinearity (Fig. 10a). This is plausible because localized shear around the splitting crack is not expected to contribute significantly to the global stiffness of the specimen. However, the decrease in specimen stiffness beyond the proportional limit is poorly replicated if shear nonlinearity is not modeled. This may be due to transverse cracks not being reproduced faithfully, because it has been established previously that the gradual decrease of specimen stiffness is due largely to transverse cracking. The splitting crack growth curve obtained without shear nonlinearity is compared in Fig. 10(b) with the curve obtained with the calibrated shear nonlinearity and the experimental curve. With shear nonlinearity the splitting crack grows even faster at relatively large applied stress, which may be surprising because shear nonlinearity should relieve strain energy and lower the crack driving force for the splitting crack. However, the greater effect is the reduction of delamination when shear nonlinearity is included, which keeps more energy available for splitting crack growth. This has also been noted recently by Van de Meer et al. (2010).

3.5.2. Influence of cohesive strengths The influences of tensile and shear cohesive strength on the specimen stiffness and splitting crack growth are shown in ^ n ¼ 50–150 MPa and Fig. 11. Fig. 11(a) shows the stress–displacement curves obtained with tensile strength varying from s shear strength from t^ s ¼ 75 to 100 MPa. ^ n . For example, when s ^ n ¼ 50 MPa, the predicted strength is 235 MPa. As s ^n The ultimate strength increases with s increases to 150 MPa, the strength is increased to 306 MPa. A smaller tensile cohesive strength results in a larger transverse crack density, and earlier development of transverse cracking near the notch tip (Fig. 11b), which shifts more load from 901-plies to the 01-plies. A larger stress concentration in the 01-plies leads to earlier fiber rupture. From the ^ n is 720 MPa. discussion presented above on calibration procedures, a reasonable estimate of the uncertainty in s The shear cohesive strength does not affect the stress–displacement relation significantly although a slight increase of stiffness is associated with shear strength increase from 75 to 85 MPa. However, the splitting crack growth rate is very sensitive to the shear strength. This is shown in Fig. 11(b): the larger the shear cohesive strength, the slower the splitting crack growth. This is attributed to the relatively small hardening modulus in the shear stress–strain relation (Fig. 6): a small increase of shear strength requires a fairly large increase in the shear strain for criticality in the elements that are to generate the splitting crack. This in turn requires a larger applied stress to propagate the splitting crack.

3.5.3. Influence of mode I fracture toughness value The mode I and mode II toughness values were varied by 50% from the baseline calibration values, with the constraint that the toughness values in modes II and III remained equal. Changing mode I toughness (GIC) around the calibration value by 750% does not significantly change the stress–strain relation or the splitting crack growth rate and has a mild influence on ultimate strength (Fig. 12). In regard to transverse ply cracking, changing GIC alone is equivalent, for the subject lay-up, to changing the specimen thickness while keeping other material parameters fixed (Eq. (12)); the effect is

f = 3.5% GIc = 250 J/m2 GIIc = 900 J/m2

300

Experimental curve (Hallet & Wisnom 2006)

200

ˆ = 150; ˆ = 75

100

ˆ = 120; ˆ = 80 ˆ = 120; ˆ = 75 ˆ = 120; ˆ = 85 ˆ = 50; ˆ = 75

0 0

0.4 0.2 Displacement, mm

0.6

12 Splitting crack length (mm)

Nominal Stress, MPa

400

f = 3.5%

ˆ = 120; ˆ = 85 (base-line)

GIc = 250 J/m2 GIIc = 900 J/m2

ˆ = 150; ˆ = 75

8

Exp.

ˆ = 120; ˆ = 75

4

ˆ = 50; ˆ = 75

ˆ = 120; ˆ = 100 ˆ = 120; ˆ = 85

0 0

100 200 Nominal Stress (MPa)

Fig. 11. Influence of cohesive strength on (a) nominal stress–displacement relation; and (b) on splitting crack growth.

300

1368

X.J. Fang et al. / J. Mech. Phys. Solids 59 (2011) 1355–1373

400

12 1.5 GIc

base-line

ˆ

1.5 GIIc

Splitting crack length (mm)

Nominal Stress, MPa

f = 3.5% ˆ

300

0.5 GIc

200 0.5 GIIc

100 Experimental curve (Hallet & Wisnom 2006b)

0

f = 3.5%

base-line

ˆ 1.5GIc

ˆ

0.5GIc

8

1.5GIIc 0.5 GIIc

4

Experimental curve (Hallet & Wisnom 2006b)

0 0

0.2

0.4 Displacement, mm

0.6

0.8

0

100 200 Nominal Stress (Mps)

300

Fig. 12. Influence of fracture toughnesses, varied around baseline values (inset), on (a) nominal stress–displacement relation and (b) on splitting crack growth.

Split tip

Split tip

Delamination front

Delaminationfront Delamination front

Split tip

Delaminationfront Delamination front

Split tip

Delaminationfront Delamination front

Fig. 13. Coupled evolution of splitting-delamination in a specimen with GIIC being reduced by 50% (0.5GIIC). Compared to the base-line results in Fig. 7, both the splitting crack growth rate and delamination rate are much higher at similar stress levels. The delamination zone is no longer a sharp triangular wedge (a, b). Rather, it propagates quickly across the specimen width bounded by the splitting crack (b, c), and quickly establishes a steady-state delamination front that propagates towards the loading edge (d). The final stage is numerically unstable under displacement controlled loading. (a) s ¼ 80 Mpa, (b) s ¼200 Mpa, (c) s ¼ 230 Mpa and (d) s ¼ 220 Mpa.

similar to changing the cohesive strength. However, the influence on ultimate strength of changing GIC by 750% is much smaller than changing the normal cohesive strength by the same proportion (Eq. (12)).

3.5.4. Influence of mode II fracture toughness value The mode II fracture toughness affects both the splitting crack, which is approximately a shear (mode II) crack, and delamination. The splitting crack propagates at a rate that is negatively correlated with the value of GIIC (Fig. 12b), as also noted by Wisnom and Chang (2000) and Hallett and Wisnom (2006a) in their study that employed pre-planted cohesive zone elements for the splitting crack and a single transverse crack. In the present study, lowering GIIC by 50% not only caused more rapid splitting crack growth (Fig. 12b), but also led to a transition in the failure sequence because of its effect on delamination. Early splitting crack growth was accompanied by rapid spreading of a delamination towards the center of the specimen (bottom edge in Fig. 4b) as well as along the split. The rapid delamination deters the development of transverse cracking, which remain few and concentrated near the notch tip with a density of  0.8 mm  1 (slightly smaller than  1.0 mm  1 for the baseline case). Shortly after the delamination spreads across the specimen, a steady-state delamination front is formed, which propagates towards the loading edge in a numerically unstable fashion. The above process is illustrated by the four snap shots in Fig. 13. In these figures, the fiberdirection strains (e11) in the 01-ply at four different stress levels are shown as contour plots. The delamination fronts and splitting crack tips are indicated. Throughout this coupled splitting and delamination process, the maximum fiberdirection strain within the delamination zones is 2.2%, well below the fiber-failure strain of 3.5% (Table 1).

X.J. Fang et al. / J. Mech. Phys. Solids 59 (2011) 1355–1373

1369

This coupled damage evolution is different from all other cases and results in a feature-rich, non-monotonic load– displacement curve (Fig. 12(a)). The initial softening of the curve is similar to other cases. However, at a stress level of about 230 MPa, the load decreases as the delamination front is established. The snap-back in the load–displacement curve corresponds to the delamination front passing the displacement recording points (dots in Fig. 4). As the delamination front and splitting crack propagate further away from the notch tip (with the splitting always well ahead of the delamination front, see Fig. 13), the section of the 01-ply that passes between the notches is increasingly decoupled from the rest of the composite and carries an increasingly complete share of the load. The stress concentration from the notch-tips disappears. While the simulation fails due to numerical instability, the anticipated specimen response and ultimate failure can be estimated as that of a single 01-ply with a width of 0.5 W (10 mm), which is indicated by the dashed line in Fig. 12, terminating at the fiber failure strain ef ¼3.5% (sf ¼ 385 MPa). 4. Discussion 4.1. Mesh requirements Mesh independence is achieved provided the element size is less than the length of the nonlinear fracture process zones of the different crack systems. When cohesive parameters are varied so that the process zone length becomes smaller, mesh refinement may be necessary to maintain accuracy. One example of this requirement was observed in the study case where GIIC was reduced by 50% from the best-calibration value: the mesh size also had to be reduced to capture the transition in failure mode and expected reduction of ultimate strength correctly. 4.2. Calibration The calibration procedures used in this paper are not ideal, but are a pragmatic use of available data for the subject ^ n , by E-glass/913 composite materials, supplemented, in the absence of experiments for calibrating the mode I parameter, s the micromechanical model of Eq. (12). The search for tests for calibrating cohesive laws, with at least the peak traction determined as well as the fracture toughness values, remains an active research topic. In some applications, it may also be useful to determine characteristics of the shape of the cohesive law; what details of the shape influence engineering fracture behavior is not yet known, but can be addressed in due course by high-fidelity simulations such as those presented here. 4.3. Statistical variance in model parameters and covariance with engineering properties By considering either the experiments that yielded calibration data or the assumptions that underlie Eq. (12), estimates have been made of the uncertainty in the calibration values of the parameters in the cohesive laws for mode I and mode II fracture. These are shown in Table 2. Since they arise from deviations in the data from tests on typical specimens, the variances are likely to be material variations, rather than errors in analysis that could be reduced by more accurate procedures. If so, they imply unavoidable variance in predicted engineering properties, which will be manifested as deviance in the test data used to qualify materials. Such deviance is essential information for design. The expected variance in predicted properties can be estimated by varying the cohesive parameters used in simulations and re-computing predictions, thus evaluating the partial derivatives

y0 @j j0 @y

ð13Þ

@j/@y, where y denotes a cohesive law parameter (GIC, etc.), y0 its best-choice calibration value, j a property (ultimate strength, etc.), and j0 its value when all cohesive law parameters take their best-choice calibration values. The partial derivatives of Eq. (13) are measures of the sensitivity of engineering properties to the material properties of fracture mechanisms. Normalization with respect to y0 and j0 allows the sensitivity measures to be independent of units. A matrix Table 2 Estimates of uncertainty in cohesive law parameters and the covariance matrix between the set of parameters and three key metrics of composite failure. Cohesive parameter

GIC GIIC s^ n t^ s ð ¼ t^ t Þ a

Sensitivity matrix ðy0 =f0 Þð@f=@yÞ, j ¼ engineering property, y ¼cohesive parameter

Estimated uncertainty in cohesive parameter (%)

10 10 15 10

Very nonlinear dependence.

Ultimate strength

Global stiffness in nonlinear regime

Splitting crack growth rate

0.1 0.25a 0.3 0.3

Small Small 0.1 0.2

0  1a 0 1

1370

X.J. Fang et al. / J. Mech. Phys. Solids 59 (2011) 1355–1373

of their estimated values appears in Table 2 for the following engineering properties: the composite ultimate strength, the global stiffness averaged over the nonlinear regime, and the growth rate of the splitting cracks averaged as the crack length grows from zero to its maximum. Other engineering propertes might be chosen, but these three cover the main features of the global performance of this material in the notched tension test. The sensitivity matrix is only qualitatively indicative of trends in performance and does not substitute for complete analysis of different cases. For variations in GIIC, in particular, the transition in failure sequence described above leads to highly nonlinear dependence in the selected engineering properties; for example, the ultimate strength is almost independent of GIIC for increasing GIIC above the best-calibration value, but drops significnatly for decreasing GIIC, when it is measured by the peak load prior to unstable behavior (Fig. 12a). The predicted variance of any quantity in experimental tests should be predicted by the product of the estimated variance in a cohesive parameter and the pertinent element of the sensitivity matrix. Combining the effects of variance in all four cohesive parameters, Table 2 predicts, for example, that the variance in composite ultimate strength should be approximately 10%. The actual variance reported for double notch tension tests of the subject composites is close to this: the strengths reported by Hallett and Wisnom (2006b) range from 243 to 291 MPa. 4.4. Transverse and delamination crack interactions In the present simulations, it is implicitly assumed that, once a transverse crack is initiated, it immediately propagates as a steady-state crack, because the in-situ strength is derived from steady-state crack analysis. For polymer matrix composites, which typically have mode II toughness several times larger than mode I toughness, this is justified. A recent 3D analysis of an individual transverse crack, which includes the influence of possible accompanying delamination cracks has shown that, if GIIC/GIC 41, the initiation stress is only 5% less than the steady-state propagation stress (Zhou et al., 2010 and works cited therein). The same analysis (Zhou et al., 2010) also showed that, as long as GIIC/GIC 41, the in-plane stress required for any delamination initiated by a transverse crack to propagate away from the transverse crack is significantly higher than the transverse crack propagation stress. Consistently, for the best-calibration cohesive parameters, the present simulations show that no delaminations were triggered by transverse cracks. However, when a GIIC was reduced by 50% in the present simulations, large scale delamination occurred with minimal transverse ply cracking. 5. Conclusions In this paper, the elemental formulation of A-FEM has been presented as a standard element that is fully compatible with any FE programs including commercial codes without access to source code. The capability of the A-FEM to account for arbitrary, nonlinearly coupled, multiple cracking systems typical in laminated composites with stress concentrators, has been demonstrated by replicating the damage systems observed in tension tests of an orthogonal double-notched tension specimen. The A-FEM formulation represents three different crack systems and fine-scale distributed shear deformation, all of which interact. Mesh independence is achieved provided the element size is less than the length of the nonlinear fracture process zones of the different crack systems. The required mesh size can be estimated in advance using analytical results for the process zone lengths. A combination of independent data and a micromechanical model were used to calibrate parameters in cohesive models that are used to represent nonlinear fracture processes and the distributed nonlinear shear deformation. With this calibration, quantitative agreement is obtained to within experimental variance between measured and predicted fracture test behavior, including the entire nonlinear stress–strain curve, the ultimate strength, the growth rates of delamination and splitting cracks, and the density of transverse ply cracks. Analysis shows that the variance in fracture properties expected for the estimated variance in material parameters matches the variance in the measured ultimate strength. Thus future elaboration of the simulations developed here can be the basis for statistical predictions of material performance. Comparative simulations show unequivocally the necessity of including the all the major damage processes explicitly to reproduce all observed damage characteristics. If incomplete nonlinear processes were included, a simulation may still be able to reproduce some aspects of the experiments, but not all in a consistent way. For example, if the shear nonlinearity is not considered, one may still obtained the stress–displacement curve and splitting crack growth curve more or less right, but the delamination is predicted incorrectly. If multiple transverse ply cracking is not included, the predicted delamination and splitting crack growth may still be close to experimental measurements, but the gradual stiffness reduction in the stress–displacement curve would be missed.

Acknowledgments The authors gratefully acknowledge the support of NASA Langley (Contract no. NNL08AA19C) and the National Hypersonic Science Center for Materials and Structures (Contract no. FA9550-09-1-0477). We thank Dr. Stephen R. Hallet for kindly providing the recalibrated stress–displacement curve in Figs. 8 and 10–12.

X.J. Fang et al. / J. Mech. Phys. Solids 59 (2011) 1355–1373

1371

Appendix A. Mixed-mode cohesive law The CZM employed in the A-FEM formulation in this study uses independent cohesive laws for the opening mode (mode I) and shearing modes (modes II and III), so that the toughnesses (total areas underneath the cohesive laws) and cohesive strengths for the three modes can be different. This type of cohesive laws as first proposed by Yang, Thouless et al. (1999, 2001) and Yang and Thouless (2001) in 2D and later extended to 3D by Yang and Cox (2005) into 3D for arbitrary interface delamination problems. The CZM accommodates separation of the total energy absorbed during fracture, Gtot, into the opening and shear components, GI, GII, and GIII: Gtot ¼ GI þ GII þ GIII

ðA:1Þ

where the separate components can be calculated by integration of the modes I–III traction–separation curves (shaded areas in Fig. A1): GI ðdn Þ ¼

Z

dn 0

GII ðds Þ ¼

Z

dt 0

GIII ðdt Þ ¼

Z 0

ds

sðd0n Þ dd0n

ðA:2aÞ

ts ðd0s Þ dd0s

ðA:2bÞ

tt ðd0t Þ dd0t

ðA:2cÞ

where dn, ds, and dt denote the mode I, II, and III displacement jumps across the interface and s, ts, and tt the corresponding mode I, II, and III cohesive tractions, respectively. Note that the energy release rate components are not independent parameters; they evolve together as a natural result of the interplay between the deformation of the bonded structural elements and the physics embedded in the traction–separation laws. A criterion is required to specify conditions for complete failure. A frequently used failure criterion is the following linear one (Wang and Suo, 1990B; Hutchinson and Suo, 1992): D ¼ GI ðdn Þ=GIC þ GII ðds Þ=GIIC þGIII ðdt Þ=GIIIC ¼ 1

ðA:3Þ

where GIC ¼ GI ðdnc Þ; GIIC ¼ GII ðdsc Þ and GIIIC ¼ GIII ðdtc Þare the mode-I, -II and -III fracture toughness in the LEFM sense. D is a damage measure ranging from 0 (undamaged) to 1 (fully debonded). When D ¼ 1, all tractions are set to zero, even though individual displacement components will not generally have reached the limiting values, dnc, dsc, or dtc, shown in Fig. A.1. In this paper we assume the same cohesive traction law for mode II and mode III, in the absence of definite experimental information about mode III delamination fracture in the literature. Finally, since in this mixed-mode cohesive model, the normal and tangential tractions are not coupled, the tangential stiffness matrix can be obtained as 2 3 kn ðdn Þ 0 0 6 0 0 0 7 ks ðds Þ Dcoh ¼ 4 ðA:4Þ 5 0 0 kt ðdt Þ where kn ðdn Þ, ks ðds Þ, and kt ðdt Þ are the instantaneous slopes in the corresponding mode I–III traction–separation laws.

n ˆ n ˆ2 

Mode -I mode-I GI n1

-sc -s2

n2 s ˆs ˆs2 -s1

nc

n t

Mode -II mode-II GII s1 -ˆs2 -s

-tc s2

sc s

-t2

ˆt ˆt2 -t1

Mode -III GIII III t1 -ˆ t2 -ˆ t

t2

tc t

Fig. A1. The modified mode-dependent cohesive zone model from Yang and Cox (2005).

1372

X.J. Fang et al. / J. Mech. Phys. Solids 59 (2011) 1355–1373

References Belytschko, T., Mo¨es, N., et al., 2001. Arbitary discontinuities in finite elements. International Journal for Numerical Methods in Engineering 50, 993–1013. Camanho, P.P., Davila, C.G., et al., 2006. Prediction of in situ strengths and matrix cracking in composites under transverse tension and in-plane shear. Composite Part A: Applied Science and Manufacturing 37, 165–176. Carpinteri, A., Ferro, G., 2003. Fracture assessment in concrete structures. In: Milne, I., Ritchie, R.O., Karihaloo, B. (Eds.), Concrete Structure Integrity. Elsevier Science, Amsterdam, pp. 7.  R., 1996. Bridged versus cohesive crack in the flexural behavior of brittle matrix composites. International Journal of Fracture 81, Carpinteri, A., Massabo, 125–145. Case, S.W., Reifsnider, K.L., 1999. MRLife 12 Theory Manual—Composite Materials. Materials Response Group, Virginia Polytechnical Institute and State University. Chessa, J., Wang, H., et al., 2003. On the construction of blending elements for local partitioning if unity enriched finite elements. International Journal for Numerical Methods in Engineering 57, 1015–1038. Choi, H.Y., Chang, F.K., 1992. A model for predicting damage in graphite/epoxy laminated composites resulting from low-velocity point impact. Journal of Composite Materials 26, 2134–2169. Cox, B.N., Dadkhah, M.S., et al., 1994. Failure mechanisms of 3D woven composites in tension, compression, and bending. Acta Metallurgica et Materialia 42, 3967–3984. Cox, B.N., Marshall, D.B., 1994. Concepts for bridged cracks in fracture and fatigue. Acta Metallurgica et Materialia 42 (2), 341–363. Cox, B.N., Spearing, S.M., et al., 2008. Practical challenges in formulating virtual tests for structural composites. In: Camanho, P.P., Da´vila, C.G., Pinho, S.T., Remmers, J.J.C. (Eds.), Mechanical Response of Composites, 10. Springer Science and Business Media, Dordrecht, pp. 57–75. Cox, B.N., Yang, Q.D., 2006. In quest of virtual tests for structural composites. Science 314, 1102–1107. Cui, W.C., Wisnom, M.R., 1992. Contact finite-element analysis of 3-point and 4-point short beam bending of uniderictional composites. Composite Science and Technology 45, 323–334. Davila, C.G., Camanho, P.P., et al., 2005. Failure criteria for FPR laminates. Journal of Composite Materials 39, 323–345. de Borst, R., Remmers, J.J.C., et al., 2006. Mesh-independent discrete numerical representations of cohesive-zone models. Engineering Fracture Mechanics 73 (2), 160–177. ¨ N., et al., 2001. An extended finite element method for modeling crack growth with frictional contact. Computer Methods in Applied Dolbow, J., Moes, Mechanics and Engineering 190, 6825–6846. Dvorak, G.J., Laws, N., 1987. Analysis of progressive matrix cracking in composite laminates. II. First ply failure. Journal of Composite Materials 21, 309–329. Fries, T.P., Belytschko, T., 2006. The intrinsic XFEM: a method for arbitryr discontinuities without additional unknowns. International Hournal for Numerical Methods in Engineering 68, 1358–1385. Galessgen, E.H., Riddel, W.T., et al., 2002. Nodal constraint, shear deformation and continuity effects related to the modeling of debonding of laminates using shell elements. Computer Modeling in Engineering and Science 3 (1), 103–116. Gonzalez, C., LLorca, J., 2006. Multiscale modeling of fracture in fiber-reinforced composites. Acta Materialia 54, 4171–4181. Hallett, S.R., Wisnom, M.R., 2006a. Numerical investigation of progressive damage and the effect of layup in notched tensile tests. Journal of Composites Materials 40, 1229–1245. Hallett, S.R., Wisnom, M.R., 2006b. Experimental investigation of progressive damage and the effect of layup in notched tensile tests. Journal of Composite Materials 40, 119–141. Hansbo, A., Hansbo, P., 2002. An unfitted finite element method, based on Nitsche0 s method, for elliptic interface problems. Computational Methods in Applied Mechanics and Engineering 191, 5537–5552. Hansbo, A., Hansbo, P., 2004. A finite element method for the simulation of strong and weak discontinuities in solid mechanics. Computational Methods in Applied Mechanics and Engineering 193, 3523–3540. Hutchinson, J.W., Suo, Z., 1992. Mixed mode cracking in layered materials. Advances in Applied Mechanics 29, 63–191. Huynh, D.B.P., Belytschko, T., 2009. The extended finite element method for fracture in composite materials. International Journal for Numerical Methods in Engineering 77 (214–239), 2009. Iarve, E.V., Mollenhauer, D., et al., 2005. Theoretical and experimental investigation of stress redistribution in open-hole composite laminates due to damage accumulation. Composites: Part A 36, 163–171. Jelf, P.M., Fleck, N.A., 1994. The failure of composite tubes due to combined compression and torsion. Journal of Materials Science 29, 3080. Kumosa, M., Odegard, G., et al., 2002. Comparison of the 7 45 tensile and Iosipescu shear tests for woven fabric composites. Journal of Composites Technology and Research 24 (1), 3–15. Lapczyk, I., Hurtado, J., 2007. Progressive damage modeling in fiber-reinforced materials. Composites Part A 38, 2333–2341. Laurin, F., Carrere, N., et al., 2007. A multi-scale progressive failure approach for composite laminates based on thermodynamical viscoelastic and damage models. Composites Part A 38, 198–209. Ling, D.S., Yang, Q.D., et al., 2009. An augmented finite element method for modeling arbitrary discontinuities in composite materials. International Journal of Fracture 156, 53–73. Maimi, P., Camanho, P.P., et al., 2007. A continuum damage model for composite laminates: Part I—constitutive model. Mechanics of Materials 39, 897–908. Matzenmiller, A., Lubliner, J., et al., 1995. A constitutive model for anisotropic damage in fiber-composites. Mechanics of Materials 20, 125–152. Mergheim, J., Kuhl, E., et al., 2005. A finite element method for the computational modeling of cohesive cracks. International Journal for Numerical Methods in Engineering 63, 276–289. Mergheim, J., Kuhl, E., et al., 2007. Towards the algorithmic treatment of 3D strong discontinuities. Communications in Numerical Methods in Engineering 23, 97–108. ¨ N., Belytschko, T., 2002. Extended finite element method for cohesive crack growth. Engineering Fracture Mechanics 69, 813–833. Moes, Pagano, N.J., Schoeppner, G.A., 2000. Delamination of Ploymer Matrix Composites: Problems and Assessment. Comprehensive Composite Materials, 2. A. Z. C. Kelly, Elsevier Science, Oxford, pp. 433–528. Pinho, S.T., Robinson, P., et al., 2006. Fracture toughness of the tensile and compressive fibre failure modes in laminated composites. Composite Science and Technology 66, 2069–2079. Remmers, J.J.C., de Borst, R., et al., 2003. A cohesive segments method for the simulation of crack growth. Computational Mechanics 31 (1–2), 69–77. Shahwan, K.W., Waas, A.M., 1997. Non-self-similar decohesion along a finite interface of unilaterally constrained delaminations. Proceedings of the Royal Society of London A 453, 515–550. Song, J.H., Areias, P.M.A., et al., 2006. A method for dynamic crack and shear band propagation with phantom nodes. International Journal of Numerical Methods in Engineering 67, 868–893. Spearing, S.M., 2010. Private communication. Talreja, R., 2006. Mulitscale modeling in damage mechanics of composite materials. Journal of Material Science 41, 6800–6812. Tay, T.-E., 2003. Characterization and analysis of delamination fracture in composites: an overview of developments from 1990 to 2001. Applied Mechanics Reviews 56 (1), 1–32. Thouless, M.D., 1990. Crack spacing in brittle films on elastic substrates. Journal of American Ceramic Society 73, 2144–2146. Turon, A., Camanho, P.P., et al., 2006. A damage model for the simulation of delamination in advanced composites under variable-mode loading. Mechanics of Materials 38, 1072–1089.

X.J. Fang et al. / J. Mech. Phys. Solids 59 (2011) 1355–1373

1373

Van de Meer, F.P., Oliver, C., et al., 2010. Computational analysis of progressive failure in a notched laminate including shear nonlinearity and fiber failure. Composite Science and Technology 70, 692–700. Van de Meer, F.P., Sluys, L.J., 2009a. Continuum models for the analysis of progressive failure in composite laminates. Journal of Composite Materials 43, 2131–2156. Van de Meer, F.P., Sluys, L.J., 2009b. A phantom node formulation with mixed mode cohesive law for splitting in laminates. International Journal of Fracture 158, 107–124. Wang, J.S., Suo, Z., 1990. Experimental determination of interfacial toughness using brazil-nut-sandwich. Acta Metallurgica 38, 1279–1290. Wisnom, M.R., 1994. The effect of fibre rotation in 7 451 tension tests on measured shear properties. Composites 26, 25–32. Wisnom, M.R., Chang, F.-K., 2000. Modelling of splitting and delamination in notched cross-ply laminates. Composites Science and Technology 60, 2849–2856. Wisnom, M.R., Jones, M.I., 1996. Size effects in interlaminar tensile and shear strength of unidirectional glass fibre/epoxy. Journal of Reinforced Plastics and Composites 15, 2–15. Xie, D., Amit, G., et al., 2006. Discrete cohesive zone model to simulate static fracture in 2D tri-axially briaded carbon fiber composites. Journal of Composite Materials 40, 2025–2046. Yang, Q.D., Cox, B.N., 2005. Cohesive zone models for damage evolution in laminated composites. International Journal of Fracture 133, 107–137. Yang, Q.D., Cox, B.N., et al., 2011. Virtual testing for advanced aerospace composites: advances and future needs. Journal of Engineering Materials and Technology, 133 011002, 1–6. Yang, Q.D., Thouless, M.D., 2001. Mixed mode fracture of plastically-deforming adhesive joints. International Journal of Fracture 110, 175–187. Yang, Q.D., Thouless, M.D., et al., 1999. Numerical Simulations of adhesively-bonded beams failing with extensive plastic deformation. Journal of the Mechanics and Physics of Solids 47, 1337–1353. Yang, Q.D., Thouless, M.D., et al., 2001. Elastic–plastic mode-II fracture of adhesive joints. International Journal of Solids and Structures 38, 3251–3262. Zhou, Z.Q., Fang, X.J., et al., 2010. The evolution of a transverse intra-ply crack coupled to delamination cracks. International Journal of Fracture 165, 77–92. Zienkiewicz, O.C., Taylor, R.L., 2005. The Finite Element Method for Solid and Structural Mechanics. Elsevier Butterworth-Heinemann, Oxford.