Numerical characterisation of the fibre–matrix interface crack growth in composites under transverse compression

Engineering Fracture Mechanics 75 (2008) 4085–4103

Contents lists available at ScienceDirect

Engineering Fracture Mechanics journal homepage: www.elsevier.com/locate/engfracmech

Numerical characterisation of the fibre–matrix interface crack growth in composites under transverse compression E. Correa *, V. Manticˇ, F. París Group of Elasticity and Strength of Materials, School of Engineering, University of Seville, Camino de los Descubrimientos s/n, 41092 Seville, Spain

a r t i c l e

i n f o

Article history: Received 7 May 2007 Received in revised form 5 March 2008 Accepted 10 March 2008 Available online 14 March 2008 Keywords: Composites Boundary element analysis Micromechanics Crack growth Interface fracture

a b s t r a c t Inter-fibre failure under compression transverse to the fibres is studied at micromechanical level. Interfacial fracture mechanics concepts, associated to both the open model and the contact model, are applied. A numerical study is performed using the boundary element method aimed at explaining the origin and evolution of the damage at micromechanical level, considered as fibre–matrix interface cracks. Assuming that the damage starts as small debonds originated by shear stresses at the position where their maximum values are reached, it has been found that the crack shows different morphologies at both tips: an open one and a closed one with a large contact zone. Then the interface crack grows unstably in mixed mode only on the open tip side until this growth changes to stable, once the crack closes at this tip, with the generation of a contact zone. Ó 2008 Elsevier Ltd. All rights reserved.

1. Introduction Knowledge of damage mechanisms at micromechanical level in composite materials is essential for the elaboration of physically based failure criteria to be applied in the design of composite structures [1]. Fibre reinforced composite materials are usually designed to work mainly in the direction of the fibres (longitudinal direction) so that they are particularly prone to failure in the transverse direction. In particular, the transverse failure stands out in the case of multidirectional laminates containing 90° plies and in laminates subjected to impact loads. With reference to the first situation, the inclusion of 90° plies in multidirectional laminates is very common due to the extra stiffness they confer on the whole structure and their capacity for preventing it from splitting. In return, they are the first plies within the laminates to show cracks due to the predominant load transverse to the fibres of these plies, the development of their failure then becoming very significant in the global damage mechanism of the whole laminate. With reference to the second situation, impact loads produce stresses that are transferred along many directions within the laminate, the transverse failure then being possible for particular plies. The mechanism of failure associated to transverse loads is commonly known as inter-fibre failure (also known as matrix failure). In previous works by París and co-workers [2–5], it was demonstrated that, under the particular case of tension loads, this mechanism of failure is characterised by the appearance of small debonds (interface cracks) at the fibre–matrix interfaces that grow along the interfaces. After these interface cracks have grown to a certain length at the interfaces, they change their direction of propagation, kinking into the matrix and propagating unstably through it, thus leading to the macro failure of the ply. The conclusions of these works lead the interface crack growth to be considered as an essential step in the whole damage process, its role being even more emphasised in the case of cyclic loads [6].

* Corresponding author. Tel.: +34 954487299; fax: +34 954461637. E-mail addresses: [email protected], [email protected] (E. Correa). 0013-7944/$ - see front matter Ó 2008 Elsevier Ltd. All rights reserved. doi:10.1016/j.engfracmech.2008.03.005

4086

E. Correa et al. / Engineering Fracture Mechanics 75 (2008) 4085–4103

Nomenclature b Dundurs bimaterial parameter l shear modulus of the material E Young modulus of the material m Poisson coefficient of the material K(l) = K1 + iK2 stress intensity factor (SIF) wK(l) SIF phase angle l characteristic length e oscillatory index of the interface crack G(Da) = GI + GII energy release rate (ERR) Da virtual crack extension wG(Da) energetic phase angle Gc(w) critical ERR k fracture mode sensitivity parameter initial length of the debond ai

In the present work, the study of inter-fibre failure under compression transverse to the fibres is undertaken, focusing on the analysis of debonds on fibre–matrix interfaces, assuming that their initiation and growth constitute not only the first step in the macro-failure but also a basic part of the whole mechanism, as occurs for the tension case. Due to some special features, such as the curved form of the interface and the expected contact between fibre and matrix surfaces of the crack (in view of the compressive character of the applied load), interface crack propagation under transverse compression represents a complex problem whose study requires the application of both interfacial fracture mechanics models: the open model [7,8] and the contact model [9]. The structure of the present paper is described in what follows: a revision of the current state of the art of interfacial fracture mechanics is presented in Section 2, introducing all the concepts applied in the subsequent analysis of the particular problem under study. The main features of the boundary element model (involving a contact algorithm) employed for the analysis are described in Section 3. In Section 4, the initiation of interface cracks at the fibre–matrix interface is assessed, their further growth being studied in Section 5. The results provided by the numerical study are evaluated and analysed applying the energetic approach of interfacial fracture mechanics. A discussion of the relevance of the results obtained for actual composites is presented in Section 6. The conclusions obtained, presented in Section 7, allow the damage at the first initial stage of failure, the interface crack growth under transverse compression, to be characterised. 2. Background to the interfacial fracture mechanics Consider a crack located at the interface between two homogeneous isotropic linear elastic materials perfectly bonded along the rest of the interface and subjected to a plane strain state, see Fig. 1 where the neighbourhood of a locally flat crack tip is shown. At the interface crack both materials may separate or maintain contact. An elastic state at the interface crack tip is characterized by the Dundurs bimaterial parameter [10] b¼

l1 ðj2  1Þ  l2 ðj1  1Þ ; l1 ðj2 þ 1Þ þ l2 ðj1 þ 1Þ

jbj 6 0:5;

ð1Þ

where jk = 3  4mk and lk ¼ Ek =2ð1 þ mk Þ is the shear modulus, Ek being the Young modulus and mk the Poisson coefficient of material k ¼ 1; 2 (fibre or matrix in the present study).

Material 1 y r θ

x

Crack tip

Material 2 Fig. 1. Local reference system at the interface crack tip.

E. Correa et al. / Engineering Fracture Mechanics 75 (2008) 4085–4103

4087

Two basic models of interface cracks developed in the fracture mechanics framework will be described in what follows: the open model by Williams [7], assuming traction free crack faces at the crack tip, and the (frictionless) contact model by Comninou [9], assuming a contact zone adjacent to the crack tip during a load application. A recent comprehensive review of both models can be found in Manticˇ et al. [11]. 2.1. Open model Introducing a characteristic length l due to Rice [8], the near-tip traction vector ahead of the crack tip along the interface (h = 0) and the near-tip relative displacement across the crack, Dui(r) = ui(r, h = p)  ui(r, h = p) (i = x, y), can be expressed in terms of the complex stress intensity factor (SIF) KðlÞ ¼ K 1 ðlÞ þ iK 2 ðlÞ ¼ jKjeiwK ðlÞ (the argument l will frequently be omitted for the sake of brevity) for r ? 0 as rffiffiffiffiffiffi  K r ie 8 K r r ie ðryy þ irxy Þðr; h ¼ 0Þ ¼ pffiffiffiffiffiffiffiffi ; ðDuy þ iDux ÞðrÞ ¼ ; ð2Þ  l 1 þ 2ie coshðpeÞE 2p l 2pr where i is the imaginary unit, e is the so-called oscillation index of the interface crack e¼

1 1b ln ; 2p 1 þ b

ð3Þ

and E ¼ 2E01 E02 =ðE01 þ E02 Þ is the harmonic mean of the Young moduli E0k ¼ Ek =ð1  m2k Þ. An unexpected basic aspect of the near-tip elastic solution of this model in (2) for b 6¼ 0 6¼ e is that due to the presence of the term (r/l)ie stresses and displacements start to oscillate when the crack tip is approached (r ? 0). As a consequence of these displacement oscillations, an infinite number of regions where the crack faces interpenetrate is predicted by this solution. The size of the zone where these physically non-admissible interpenetrations occur is frequently very small. According to Rice [8], when this size is sufficiently small in comparison with the smallest characteristic length of the specimen (e.g. crack length or an adjacent layer thickness), then the open linear elastic model is adequate for interface crack growth predictions, such a situation being referred to hereinafter as small-scale contact (SSC) case (cf. Section 2.2). As follows from (2), the phase angle wK(l) = argK(l) can be evaluated using the traction vector components at the distance l ahead of the crack tip as wK ðlÞ ¼ argðrxy þ iryy Þðr ¼ l; h ¼ 0Þ;

ð4Þ

where arg represents the argument function of a complex number. Notice that tan wK ðlÞ ¼ K 2 ðlÞ=K 1 ðlÞ ¼ rxy =ryy ðl; 0Þ. Thus, wK(l) represents an l-dependent measure of fracture mode mixity, although usually varying only very slowly with l following 0 0 the rule wK ðl Þ ¼ wK ðlÞ þ e ln ll . The distance qi from the crack tip to the first interpenetration point can be estimated by the general expression recently introduced by Graciani et al. [12], see also Hills and Barber [13]: qi ffi l exp



   1 p  wK sgnðeÞ þ arctanð2jejÞ 2n  jej ; 2

ð5Þ

where n is the integer giving the largest value of qi lower than the crack length. The total strain energy release rate (ERR) due to a crack extension along the interface can be evaluated generalising the classical virtual crack closure technique (VCCT) [14], to the oscillating near-tip elastic field (2), giving for a small but finite crack extension Da: Z Da 1 GðDaÞ ¼ GI ðDaÞ þ GII ðDaÞ ¼ rjy ðr; 0ÞDuj ðDa  rÞdr ðj ¼ y; xÞ: ð6Þ 2Da 0 The limit G = limDa ? 0G(Da) exists and can be expressed, according to the classical relation by Malyshev and Salganik [15], by G¼

jKj2 2

cosh ðpeÞE

ð7Þ

;

in terms of the SIF modulus. The previous relation can be seen as a generalization to interface cracks of the Irwin’s relation [14] between the ERR (a global fracture parameter) and the SIF (a local fracture parameter) for a crack in homogeneous materials. However, there is no limit value of the components of ERR, GI(Da) and GII(Da), for Da ? 0 and e 6¼ 0. A simple explicit expression of their values for small Da in terms of wK was deduced by Manticˇ and París [16]:

    Da ;e ; ð8Þ GI;II ðDaÞ ¼ 0:5G 1  FðeÞ cos 2 wK þ w0 l where the plus and minus signs, respectively, refer to GI(Da) and GII(Da), and the amplitude and phase shift functions may be approximated by the following series expansions (see [16] for exact expressions of these functions): FðeÞ ¼ 1 þ ðp2 =3  2Þe2 þ Oðe4 Þ;

w0 ðDa=l; eÞ ¼ e ln Da=4el þ ðfð3Þ þ 4=3Þe3 þ Oðe5 Þ;

ð9Þ

4088

E. Correa et al. / Engineering Fracture Mechanics 75 (2008) 4085–4103

: : f(3) ¼ 1.202 being the Apéry’s constant and e ¼ 2.718 the base of the natural logarithm. As follows from (8) and (9), the energetic angle wG(Da) defined by tan2 wG ðDaÞ ¼

GII ðDaÞ GI ðDaÞ  GII ðDaÞ or equivalently by cos 2wG ðDaÞ ¼ GI ðDaÞ GI ðDaÞ þ GII ðDaÞ

ð10Þ

represents a Da-dependent measure of fracture mode mixity, although usually varying only very slowly with Da (the argument Da will frequently be omitted for the sake of brevity). The following expression, obtained from (8) and (10), can be used to evaluate wK in terms of the energetic fracture parameters [16]:

 GI  GII ; ð11Þ w0K ¼ 0:5 arccos½FðeÞ1 cosf2wG g ¼ 0:5 arccos FðeÞ1 GI þ GII where 0 6 w0K ¼ jwK þ w0 þ npj 6 p=2 with n being an integer number. Note that (11) particularizes to |wK| = wG for b = 0 = e. 2.2. Contact model The contact model of Comninou [9] tries to overcome the inconsistencies of the open model. Allowing a frictionless contact between the crack faces, a physically correct solution with one contact zone at the crack tip is obtained when b 6¼ 0. Typically, the extent of this contact zone qc is smaller than the size of the interpenetration zone qi in the open model (5). A mechanical explanation of this fact was introduced by París et al. [2], and the following estimation of qc by means of qi was deduced by Hills and Barber [13]:

 arctanð2eÞ : ð12Þ qc ffi qi 4 exp  e Due to the presence of a near-tip contact, the interface crack grows in pure Mode II. Thus, the near-tip singular elastic state is uni-parametric, being governed by one multiplicative constant represented by the SIF K CII . The near-tip traction vector and the near-tip relative displacement across the crack are expressed for r ? 0 as rffiffiffiffiffiffi K CII bK CII 8K CII r : ð13Þ rxy ðr; 0Þ ¼ pffiffiffiffiffiffiffiffi ; ryy ðr; pÞ ¼  pffiffiffiffiffiffiffiffi 6 0; Dux ðrÞ ¼ 2 2pr 2pr cosh ðpeÞE 2p An important consequence of the inequality in (13), implied by a requirement of near-tip compressive stresses between crack faces, is that bK CII P 0, the sign of K CII depending only on b and being independent of the far-field load direction. Hence, the only allowed direction of the near-tip relative slip is defined by the relation bDux(r) P 0 for r ? 0. When the global imposed shear loading agrees with this intrinsically allowed slip direction a relatively large near-tip contact zone may take place. However, when the applied global load tends to originate slip opposite to the allowed near-tip slip direction, only a very small contact zone, typically of subatomic size, is present in the analytic solution at this tip (sliding in the locally allowed direction) with an adjacent crack opening, sometimes referred to as a ‘bubble’ (with the crack faces sliding in the direction imposed by the load). This behaviour of an interface crack solution was first studied by Comninou and Schmueser [17] and Gautesen and Dundurs [18] and further discussed by Leguillon [19] and Audoly [20]. According to these studies, the contact zones at both crack tips are separated by an intermediate zone of a crack opening. The effect of shear is that a large contact zone is produced at one crack tip and an extremely small one at the other crack tip, Fig. 2. As a consequence, the asymptotic singular solution of the contact model in presence of the ‘bubble’ extremely close to the crack tip becomes physically meaningless, such situations definitely fulfill SSC conditions, no near-tip contact zone is observable in experiments and, thus, the locally open model is the appropriate one (in spite of its inconsistencies) for the analysis and prediction of crack behaviour in such situations. Evaluating the VCCT integral in (6) using the near-tip solution from (13), taking j = x only, yields a relation analogous to (7) where G is replaced by GII and |K| by K CII . 2.3. Interface crack propagation The propagation of an interface crack by its further extension along the interface may happen either under SSC conditions in mixed fracture mode or, in presence of a physically relevant near-tip contact zone, in pure shear mode. Then, the standard criterion for the interface crack growth can be formulated on the energetic basis comparing the corresponding ERR for a load level, G, with the critical ERR, Gc [21], as expressed in the following inequality: G P Gc :

ð14Þ

In the SSC case, a strong dependence on the mode mixity of Gc(wK) has been observed in extensive experiments [22–24], wK (or wG) then being an important parameter governing interface crack growth. The following phenomenological law for Gc(wK) [21], is considered to be representative of a large number of bimaterial systems: Gc ðwK Þ ¼ G1c ð1 þ tan2 ð1  kÞwK Þ;

ð15Þ

E. Correa et al. / Engineering Fracture Mechanics 75 (2008) 4085–4103

4089

External load Material 1-Stiff(Fibre) Allowed near-tip slip direction

Not allowed near-tip slip direction

Large near-tip contact zone

Extremely small near-tip contact zone

“Bubble” Material 2-Compliant (Matrix)

Mat. 1

External load Local sliding direction

Mat. 2

Fig. 2. Morphology of the interface crack under shear external load.

where G1c is the Mode I critical interface ERR (also called separation energy, associated to the minimum value of Gc(wK)) and k is a fracture mode-sensitivity parameter, e.g. typical range 0.2 6 k 6 0.3 characterises interfaces with moderately strong fracture mode dependence. In the case of a physically relevant near-tip contact zone and assumed frictionless contact the critical ERR could be considered a constant due to the uni-parametric character of the near-tip solution (13). Nevertheless, in presence of a physically relevant contact zone it might be questionable to neglect the effect of frictional dissipation in this zone during crack growth, which reduces the amount of energy available for this growth. As follows from the studies by Stringfellow and Freund [25] and Sun and Qian [26], where a modification of the standard VCCT (6) was introduced to evaluate the energy available for crack growth in presence of friction, the effect of friction can be taken into account in two ways: to evaluate the reduced ERR (according to Sun and Qian [26]) and compare it with a constant intrinsic critical ERR, or to consider a critical apparent-ERR incremented by the friction dissipation term. For a not very large near-tip contact zone, a rough estimation of the critical apparent-ERR can be obtained from (15), where wK is evaluated from tractions ahead of the crack tip where compressions can be expected, which implies values of |wK| > 90°. 3. Single fibre model The numerical analysis is performed using a numerical tool based on the collocational boundary element method (BEM) [27], and developed in [28]. This BEM code permits a numerical analysis of the plane elastic problem considering contact and interface cracks to be carried out, in a similar way to that described in detail by Graciani et al. [29] for axisymmetric problems. The BEM model used in the study carried out corresponds to the problem shown in Fig. 3 and consists of a single fibre configuration that, under the hypothesis of plain strain, incorporates a circular crack situated at the fibre–matrix interface, this model being able to analyse its growth caused by the transverse compressive load applied at the external boundary of the surrounding matrix. Only stresses associated to the external load are considered in this paper. In any case, the possible presence of curing stresses was demonstrated [30] to have no qualitative influence on the shape of the evolution of the ERR with the debonding angle in fibre–matrix interface cracks under far-field tension. The BEM model permits the detection of possible contact between fibre and matrix along the debond. In order to achieve a high accuracy in the numerical results, BEM meshes strongly refined towards the crack tips are employed, the size of the smallest element located at both crack tips being 3.5  107a, where a is the fibre radius. All the boundaries in this study have been modelled with continuous linear elements [27]. To characterize the problem from the fracture mechanics point of view, the ERR will be used. Adapting the VCCT, Eq. (6), to the present circumferential crack case, when the crack propagates from a certain crack tip position given by the angle a to a + Da (Da being much smaller than the crack length), gives: Z Da 1 Gða; DaÞ ¼ frrr ða þ hÞDur ða  Da þ hÞ þ rrh ða þ hÞDuh ða  Da þ hÞgdh; ð16Þ 2Da 0 where rrr and rrh represent, respectively, radial and shear stresses along the interface, and Dur and Duh the associated relative displacements of the crack faces. In this study, the integration range chosen, Da, takes the value of 0.5°. Within this range the distribution of nodes involved in the ERR calculations is similar to that already employed by París et al. [2].

4090

E. Correa et al. / Engineering Fracture Mechanics 75 (2008) 4085–4103

3

σ0

Upper tip

α

2

a FIBRE

Lower tip MATRIX

σ0 Fig. 3. Single fibre model with an interface crack.

Table 1 Elastic properties of materials Material

Poisson coefficient

Young modulus

Matrix (epoxy) Fibre (glass)

mm = 0.33 mf = 0.22

Em = 2.79  109 Pa Ef = 7.08  1010 Pa

The bimaterial system considered in this paper corresponds to a glass–epoxy system. In a former paper [5], it was clarified that the type of phenomenon that takes place in these micromechanical configurations is much more dominated by the geometry than by the properties of the materials, which only slightly alter the results quantitatively. The elastic properties of the materials chosen are listed in Table 1. Finally, in order to completely characterize the bimaterial system considered, Dundurs bimaterial parameter b [10] takes the value: b = 0.229, while the oscillatory index e = 0.074. The geometry is defined by the radius of the fibre, a = 7.5  106 m, and the dimension of the external boundary of the matrix side has been taken equal to 103 m. Dimensionless results for ERR will be presented in all cases. These  dimensionless  m r20 ap, where values of ERR, G, are obtained, following Toya’s approach [31,32], by dividing the values of G by G0 ¼ 1þj 8lm m m m j = 3  4m , l is the shear modulus of the matrix and r0 denotes the absolute value of the applied compression. 4. Damage initiation at the interface The analysis of the stress state around a fibre with an undamaged interface has already proved to be a way to locate the initiation of damage at the interface [5]. Thus, assuming the situation of an external compression acting over an initially undamaged interface, Fig. 4, an analysis of the stress situation surrounding the fibre has been performed making use of Goodier’s analytical expressions [33], which is appropriate in view of the relative dimensions of the fibre and the surrounding matrix. The results obtained are presented in Fig. 5, where radial stresses, rrr, and shear stresses, rrh, at the interface are plotted versus the angular position, a (due to symmetry, only results for one half of the fibre are presented). These curves confirm the existence of zones of maximum shear stresses, located at a = 45°, 135°, 225° and 315°. The radial stress shows an also predictable compressive character affecting almost the entire boundary of the fibre, excepting an area of reduced size centred at a = 0° and 180°, where tensions of small value are acting. It can also be observed from the figure that the maximum shear stress is about 7.5 times higher than the maximum tensile radial stress. Although the application of the shear strength concept in engineering is quite common, it is often difficult to find good tabulations in the literature on shear strength. As a conservative relationship it is recommended [34], to take for metals the value of the shear strength as 40% of the value of the tensile strength. With reference to epoxy resin, the shear strength values are usually a little bit lower than the tensile ones [35], and thus, considering both strengths to be of the same order leads

4091

E. Correa et al. / Engineering Fracture Mechanics 75 (2008) 4085–4103

σ0

perfect bonding

α

FIBRE

MATRIX

σ0 Fig. 4. Single fibre configuration for an undamaged interface.

1.5

σ rθ /σ 0 , σ rr / σ 0

assumed length of the original debond

σ rθ srq σ rr srr

1

possible positions of the assumed debonds

0.5

0 0 -0.5

15

30

45

60

75

90

105

120

135

150

165

180

α (º)

-1

-1.5 Fig. 5. Radial stresses, rrr, and shear stresses, rrh, at an undamaged fibre–matrix interface versus the angular position, a.

to a conservative estimation. Based on this and taking into account the results presented in Fig. 5, it can be concluded that the maximum shear stress/maximum tensile stress ratio is, in any case, higher than the shear strength/tensile strength ratio. All this information makes it reasonable to assume the location of the origin of the first debonds at the interface to be at the areas of maximum shear stress, in the belief that, when this stress reaches the shear strength of the interface, a debond is produced. The length of this initial debond will depend on the shape of the shear stress curve at the neighbourhood of its maximum, so that the flatter this curve is around the peak, the longer the extension of the initial damage will be. In the present case shear stress remains quite constant in a zone corresponding to an arc length of 20° centred at the aforementioned peaks (in this interval the maximum difference in rrh is about 6%). Thus, it is quite reasonable to assume that when the shear stress reaches the shear strength value of the interface it will produce a debond of an order of between 10° and 20°, starting from which interface fracture mechanics is able to control the growth of such a defect. Based on the former reasoning and choosing arbitrarily from the four possible areas of damage initiation (a = 45°, 135°, 225°, 315°), a 10° length debond centred at 135° has been considered as the smallest initial debond whose propagation is studied in the next section. The consideration of different initial debonding angles, in the range considered, does not affect in any case the conclusions obtained in the subsequent analysis, as will be seen in Section 5.3.1.

4092

E. Correa et al. / Engineering Fracture Mechanics 75 (2008) 4085–4103

5. Numerical analysis of the fibre–matrix interface crack 5.1. Characterisation of the interface crack Based on the conclusions of the previous section, a BEM model, Fig. 3, has been developed in order to study the propagation along the interface of an initial 10° length debond centred at 135°. This propagation is analysed in terms of the ERR, calculated for both crack tips (introduced as upper tip and lower tip in Fig. 3) at all steps of growth considered in this analysis. The ERR comparison between both crack tips will make it possible to reach a decision about the probable direction of growth along the interface, which could be, in principle, either from the upper tip or from the lower tip. It is possible, from a qualitative point of view, to make a prediction of the local morphology of the crack near its tips in the light of the contact model, as explained in Section 2.2 (Fig. 2). In this sense, a physically relevant contact zone is predicted at the upper crack tip, whereas a physically meaningless contact zone is predicted at the lower crack tip. A very fine discretization is employed near both crack tips, the length (subatomic) of the smallest boundary element adjacent to the crack tip being 2.6  1012 m. Therefore, the discretized model carried out is perfectly prepared to allow the contact zone appearing at the crack tips to be always detected by the BEM contact algorithm when this zone is physically meaningful (from the continuum mechanics point of view). Thus, the BEM contact algorithm will directly show which interface crack model is suitable for application at each crack tip for crack growth assessment. The obtained numerical solution will actually confirm the above qualitative predictions, as will be seen in the course of this section. In this situation, the study starts with the comparison of G for both crack tips for the initial debond of 10° length, obtaining that Gja¼140 (i.e. ERR of the lower tip) is slightly higher than Gja¼130 (i.e. ERR of the upper tip), see Table 2. The slightly negative value of GI observed in Table 2, for almost pure fracture Mode II at a = 140°, is consistent with the theoretical results for the open model proved by Manticˇ and París [16]. It was shown there that interface crack configurations (with e 6¼ 0) working predominantly in fracture Mode I or II, respectively, may have associated a negative but small value of the ERR component II or I. This result can be easily seen from the expression of GI and GII in Eq. (8) and the fact that, according to Eq. (9), F(e) > 1. It is important to notice that although both total ERR values agree in the dominant shear fracture mode character (Mode II), checking the deformed situation of the crack faces the existence of a contact zone at the upper crack tip is detected, while a small ‘bubble’ is found at the lower crack tip. This situation is illustrated by Fig. 6, where the distribution of dur along the original interface crack is represented in polar coordinates. This representation allows a comparison between the positions of the nodes of the BEM mesh at the interface corresponding to the undeformed and deformed configurations (only radial relative displacements considered) to be carried out. The existence of the contact zone at the upper crack tip and the ‘bubble’ at the lower crack tip, in agreement with the analytical predictions of the contact model for an interface crack under the action of a shear field (see Section 2.2), can be appreciated in Fig. 6. Although the total G values that appear in Table 2 would suggest that the crack may grow in both directions, the presence of an opening zone, ‘bubble’, see Section 2.2, at the lower crack tip in opposition to the large contact zone at the upper crack tip suggests the imminent appearance of a Mode I, making the direction of growth at the lower crack tip more plausible. Based on this it is assumed that if crack propagation takes place, it will occur from the lower crack tip. Based on the same considerations as those used for the original debond the growth of larger debonds is considered to take place through the lower crack tip. Complementary studies supporting this decision will be presented later in this section. The results of the G evolution at the lower crack tip during the growth of the crack from this tip are shown in Fig. 7, where the appearance, growth and latter disappearance of Mode I in this evolution can be observed. The relevant participation of Mode II at the initial stage of crack growth is also detected as well as its decreasing character as Mode I increases. Finally, Mode II changes to increasing in the last steps of growth coinciding with the disappearance of Mode I. This change in character from mixed mode to Mode II (immediately after a = 206°) redounds in the closing of the end part of the ‘bubble’ that existed at the lower crack tip from the very beginning of growth, giving rise to a contact zone of finite size at the neighbourhood of this crack tip. Referring to the evolution of the morphology of the crack as it grows along the interface from its lower tip, the numerical results obtained show that the ‘bubble’ detected at the lower crack tip of the original debond (a = 140°, Fig. 6) increases as the crack does, progressively changing its shape until finally closing at about a = 206°. Meanwhile the faces of the crack in the neighbourhood of the upper crack tip remain in contact during the whole growth process, even enlarging the contact length at the interface, the ‘contact model’ being then of application for this tip during the range of propagation of the interface

Table 2 G, GI and GII for an interface crack between a = 130° and 140° a (°)

GI GII G

130°

140°

2.83E5 2.46E2 2.46E2

8.56E5 2.47E2 2.47E2

4093

E. Correa et al. / Engineering Fracture Mechanics 75 (2008) 4085–4103

α =130º

MATRIX

FIBRE

δ ur Undeformed position of the interface

α =140º

Fig. 6. dur situation for the original debond.

0.050 0.045

G I /G 0 Serie2 G II /G0 Serie3 G /G0 Serie4

0.040 0.035 0.030 0.025 0.020 0.015 0.010 0.005 0.000 130

140

150

160

170

180

190

200

210

220

α : position of the lower crack tip (º) Fig. 7. G of the lower crack tip of the interface crack.

crack considered. Results illustrating these facts are presented in Fig. 8, where the evolution of the radial relative displacements of the crack faces, dur, for the complete extension of the interface crack at several steps of growth are plotted. This graph shows the just-mentioned evolution of the lower tip ‘bubble’ from the initial debond (a = 140°) until its closing in the last step of growth analysed (a = 212°). The contact between the crack faces at the upper tip for all steps of crack growth considered can also be observed in this figure. Recall, examining in greater depth the deformed morphology of the interface crack, as proved by Comninou [9] and Comninou and Schmueser [17] and explained in Section 2.2, that there is always an analytically predicted contact zone at the neighbourhood of the crack tips for all cracks between dissimilar materials. When referring to the lower crack tip, the contact zone predicted may be of subatomic length (being meaningless from the continuum mechanics point of view) and become numerically undetectable in spite of the very fine BEM discretization employed. The length of the interpenetrations zone of the open model at the lower crack tip for different crack lengths can be estimated by means of the relation (5). Using the values of the phase angle wK at the lower crack tip calculated through the BEM

4094

E. Correa et al. / Engineering Fracture Mechanics 75 (2008) 4085–4103

1.2E-07

α =212º 1.0E-07

δ ur /a

8.0E-08

α =180º

6.0E-08

4.0E-08

α =172º

2.0E-08

α =140º 0.0E+00 130

140

150

160

170

180

190

200

210

220

α (º) Fig. 8. dur evolution of the interface crack.

model, a semianalytical estimation of the length of the interpenetration zone qi can then be obtained for different positions of the lower crack tip given by the angle a, see Table 3. By applying the relation (12), these values of qi can be further used to obtain estimations of the length of the contact zone qc at the lower crack tip, also shown in Table 3. In addition to these semianalytical estimations of qi and qc, Table 3 also includes the numerical results for the contact zone length detected at the lower crack tip by the BEM model, denoted as qc (BEM). The estimated values of qi and qc presented in Table 3 shows that the lengths of the interpenetration zone (in the open model) and the contact zone (in the contact model) are negligible from the continuum mechanics point of view (they being of subatomic size), for the lower crack tip angles smaller than a ffi 204°. The contact zone at the lower crack tip starts to be detectable by the present BEM model for a ffi 200°, where qc is of an order of 1011 m. Therefore, the numerical solution is shown to be able to detect a contact zone not only when it is physically meaningful, but also when it has a subatomic size. It is in any case important to observe that although both approaches (semianalytical and numerical) coincide in the order of the contact zone detected, the accuracy of the semianalytical estimation decreases as the contact zone becomes physically relevant. The explanation of this fact lies in the limited scope of the application of expressions (5) and (12), since they were deduced for straight interface cracks and considering only the first singular term of the open model (2). Moreover, the dependency of qi on wK according to (5), both being defined for the open model, clearly shows the limitation of the combination of Eqs. (5) and (12) to accurately quantify the real length of the contact zone when it reaches a relevant size, since wK is not able to correctly represent the closure of the crack tip. In particular, once wK reaches the value that indicates a pure Mode II behaviour accompanied by the appearance of a physically relevant contact zone at the lower crack tip, wK stays roughly constant for larger values of a, there thus being no correlation with the increasing length of the contact zone at this crack tip. This behaviour can be clearly appreciated in Fig. 9 (Section 5.2), where a deeper study of wK is presented. The above explained facts show that the debonding crack at the lower crack tip behaves in accordance with the open model, see Section 2.1, from a = 140° up to a ffi 200°, whereas for larger a it behaves in accordance with the contact model. When referring to the upper crack tip, as was explained before, a large contact zone is maintained, for the whole range of interface crack lengths considered, its behaviour then being in accordance with the contact model. It should be stressed that the presence of this large contact zone at the upper crack tip affects considerably the local solution at the lower crack tip. In fact, it has been numerically checked by the authors that if Toya’s analytical solution [32],

Table 3 Estimated lengths of the interpenetration, qi, and contact zones, qc, at the lower crack tip Position of the lower crack tip, a (°)

qi (semianalytical) (m)

qc (semianalytical) (m)

qc (BEM) (m)

140 142 168 180 192 200 204 206

1.34E26 2.98E25 5.11E23 1.98E20 7.83E16 3.71E11 3.27E09 1.55E08

7.36E27 1.64E25 2.81E23 1.09E20 4.30E16 2.04E11 1.80E09 8.53E09

– – – – – 1.31E11 1.25E09 1.57E08

4095

E. Correa et al. / Engineering Fracture Mechanics 75 (2008) 4085–4103

100

90º

80 60 40 20 0 140 -20

position of the crack tip, α (º) 150

160

170

180

190

200

210

220

ψ G (º) YG ψ K (ψ G ) (º) YK ψ K (σ) (º) Yk(s)

-40 -60 -80

-90º

-100 Fig. 9. wG, wK(wG) and wK(r) evolution at the lower crack tip.

which considers traction free crack faces, is used to predict the value of ERR at the lower crack tip (governed by the open model), the results obtained do not provide a sufficiently accurate approximation of the ERR distribution in comparison with the more realistic BEM results shown in Fig. 7. In particular, Toya’s solution leads to a significant overestimation of ERR for short and medium size cracks. The reason for this overestimation is that in Toya’s solution there appears a large interpenetration zone at the upper tip, corresponding to the large contact zone in the BEM model, see Figs. 6 and 8. This result indicates the need to model appropriately the solution at the upper tip, even if the interest is focused on the lower tip. 5.2. Study of the interface crack propagation To describe how the damage may grow in mixed mode it is necessary [21] to have an estimation of the critical value of G, Gc ¼ Gc =G0 , which depends on the evolution of the fracture mode mixity and, consequently, on the position of the lower crack tip. Then, first of all, knowledge of the evolution of the fracture mode mixity, characterized by the local phase angles, wK or wG, with the position of the lower crack tip is required, as was explained in Section 2. Fig. 9 represents the evolution with the angular position of the lower crack tip, a, of the following phase angles: wG (calculated from Eq. (10) and taking into consideration the sign of the shear stresses ahead of the crack tip); wK calculated at a + l from stresses ahead of the crack tip using Eq. (4) (represented as wK(r) in the figure) and wK calculated from Eq. (11) (represented as wK(wG) in the figure). For the case considered l/Da has been taken equal to 0.093245, in order to obtain a vanishing phase-shift angle w0 defined in Eq. (9)2. Selection of a different l/Da ratio would only produce a w0 translation of all points of the wK curves according to Eqs. (11) and (9)2, not altering the conclusions derived from Fig. 9. It can be observed from the figure that there is an excellent agreement between the values obtained for wG, wK(wG) and wK(r). It is relevant to highlight that wK has been defined only for the open model. Thus, when a contact zone of finite size appears (a = 206°), Mode I progressively disappears, i.e. wG = 90°, and wK(wG), obtained using Eq. (11) (based on the open model near-tip solution), no longer has strict sense. It should be mentioned that slight differences between wK(r) and wK(wG), observed for a = 212° and 214°, are due to the above-mentioned increasing size of the near tip contact zone and consequently the ceasing strict validity of the open model solution. The wK(r) evolution with the angular position of the crack tip, a, can be even better observed in Fig. 10. In this figure shear stresses are represented versus radial stresses, it then being possible to represent wK as the angle of the traction vector ahead of the crack tip at the distance l, see (4). As the interface crack grows under compression, radial stresses evolve from small negative values (for the initial steps of growth of the crack, a ffi 140°) to positive values (showing a maximum around a = 180°), and decreasing until reaching negative values again, when a large contact zone develops, a = 214°. As regards shear stresses ahead of the crack tip, it can also be observed in Fig. 10 that they start acting against the allowed near-tip slip direction, causing the ‘bubble’, and change their sign around 180°, where the ‘bubble’ starts to progressively modify its shape until it closes, leading to a contact zone of physically relevant size at the tip. Fig. 11a and b represents the distribution of radial stresses and shear stresses, respectively, ahead of the crack tip for three different angular positions (a = 140°, 164°, 198°). The aforementioned differences in the sign of stresses as the crack propagates along the interface can be clearly observed in these figures. Our attention could be attracted to the tensile character of singular radial stresses close to the crack tip for all three debonds presented. This fact could perhaps be unexpected due to the compressive far-field applied, although it is in agreement with the results shown in Figs. 9 and 10. Once the evolution of the local phase angles versus the angular position of the lower crack tip is known, the evaluation of Gc can be undertaken. In this study the evolution of the Gc considered is based on the proposal by Hutchinson and Suo [21],

4096

E. Correa et al. / Engineering Fracture Mechanics 75 (2008) 4085–4103

−σ r θ /σ 0 α =214º

8 α =206º

6 4

α =204º α =198º

2

ψΚ

0 -1

0

1

2

3

4

-2

σ rr /σ 0

5

-4 α =140º

α =172º

-6 α =164º

α =152º

-8 -10

Fig. 10. rrh versus rrr ahead of the lower crack tip for different lengths of the debonding crack defined by the lower crack tip position angle.

Eq. (15). For this analysis three different values of k (in the range of typically acceptable values) have been chosen: k = 0.2, 0.25 and 0.3. In order to relate G with Gc(wK) the value of G1c, in the absence of direct experimental data, has been adjusted to achieve, for a given far-field load, the equality G = Gc(wK) at the initial stage of growth of the crack for each value of k. Once the values of G1c have been obtained for each k considered (G1c(k = 0.2)/G0 = 0.0031, G1c(k = 0.25)/G0 = 0.0044, G1c(k = 0.3)/G0 = 0.0059),

8 6 Lower crack tip position α =140º

4

Lower crack tip position α =164º

Lower crack tip position α =198º

σrr /σ 0

2 0 130

140

150

160

170

180

190

200

210

220

α (º)

-2 -4 -6 -8 8 3

σ rθ

6

Lower crack tip position α =140º

4

Lower crack tip position α =164º

2

FIBRE

σrθ /σ 0

2 0 130 -2 -4 -6

140

150

160

170

180

190

200

210

α (º)

220

Lower crack tip position α =198º

-8 Fig. 11. Stress state ahead of the crack tip for different debonds: (a) radial stress and (b) shear stress.

4097

E. Correa et al. / Engineering Fracture Mechanics 75 (2008) 4085–4103

0.050 G G(α )/G 0 σ =σ 0 G ψ Κ (α ), λ =0.2)/ G 0 c (lambda=0.2 Gc G ψ K (α ), λ =0.25) /G 0 c (lambda=0.25 Gc G ψ K (α ), λ =0.3) /G 0 c (lambda=0.3 Gc

0.045 0.040 0.035 0.030 0.025 0.020 0.015 0.010 0.005 0.000 130

140

150

160

170

180

190

200

210

220

Position of the crack tip, α (º) Fig. 12. G and Gc for the interface crack.

the Gc, Eq. (15), has then been calculated and plotted in Fig. 12, where a comparison with the total ERR values previously obtained is presented. In this representation, the l/Da ratio, necessary for the wK calculation using Eqs. (11) and (9)2, has been taken equal to 0.75  0.093245. This value makes wK = 90° to coincide with the appearance of a contact zone of relevant size at the lower crack tip, which provides a physical meaning to the l/Da selection employed. From this figure it can be deduced that, assuming, as commented before, G(a = 140°) = Gc(a = 140°), the growth of the crack will be unstable until the lower tip reaches a position defined by a ffi 206°. The end of the unstable propagation of the interface crack coincides with the appearance of a contact zone of physically relevant size at the lower crack tip. This fact would make the growth of the interface crack beyond a ffi 206° even more difficult, due to the detrimental effect that the unavoidable presence of friction, not considered in the numerical results presented, would have on the amount of energy available for the crack to continue growing, see [6] and the pertinent comments in Section 2.3. 5.3. Justification of the initial hypotheses Finally, and now that the growth of the interface crack has been studied, it is time to present the complementary arguments that justify the initial assumptions accepted at the beginning of Section 5 related to the propagation process, i.e. low sensitivity of the results to the initial length of the interface crack and propagation of the crack taking place only at the lower tip. 5.3.1. Sensitivity to the initial length of the debond In accordance with the results of Section 4, it was decided to assume an initial defect at the fibre–matrix interface in the form of a debond of 10° length centred at the angle 135°, whose propagation as an interface crack has been studied in Sections 5.1 and 5.2. The choice of the initial length of the debond, denoted by ai in what follows, was based on the shape of the shear stress distribution at the neighbourhood of its maximum, see Fig. 5, and was located within the range 10–20°. Since the results presented so far for the propagation of the interface crack correspond to the lower bound of this range, ai = 10°, a study of the sensitivity of the results to the initial length of the interface crack is required to complete the analysis presented. To this end the evolution of the interface crack has been studied again, but now assuming ai = 20°, i.e. the upper bound of the plausible range. The analysis of the evolution of the interface crack consists of two parts: the evaluation and comparison of the ERR with its critical value, and the study of the deformed morphology of the crack near its tips. Both analyses have been performed in a similar way to that described for the case ai = 10° in Sections 5.1 and 5.2. Referring to the energetic study, its results are summed up in Fig. 13, where G and Gc for the ai = 20° case are represented together with the results for ai = 10° (already appearing in Fig. 12). It is important to note that, since both curves correspond to the same material, the Gc evolution has been calculated using the same parameters in Eq. (15) for both cases (in particular k and G1c). In this sense, as the value of G1c was chosen to fulfil the equality G = Gc(wK) at the initial position of the lower crack tip for ai = 10°, the G curve for the ai = 20° case must now be calculated not for the external load r0, but for the critical load rc that adjusts the G curve to achieve the equality G = Gc(wK) for a = 145°. The results presented in Fig. 13 first show that the propagation of an interface crack coming from a larger ai occurs at a lower level of the external load. Related to the final extension of the crack along the interface, once the propagation of the initial debond starts, it maintains an unstable character up to a position of the lower crack tip a ffi 200° for the ai = 20° case. As can be observed in Fig. 13, this position angle is slightly smaller than that obtained for the ai = 10° case (a ffi 206°).

4098

E. Correa et al. / Engineering Fracture Mechanics 75 (2008) 4085–4103

0.050 G( Gα )/G0 σ =σ0 (α i=10º) Gc(ψK((α α ), λ =0.25) /G0 (α i=10º) G( G (α α )/G60 σ =σc (αii = 20º) Serie6 Gc(ψK((7 αα ), λ =0.25) /G0 (αi =20º)

0.045 0.040 0.035 0.030 0.025 0.020 0.015 0.010 0.005 0.000 130

140

150

160

170

180

190

200

210

220

Position of the crack tip, α (º) Fig. 13. G and Gc for the interface crack in the cases ai = 10° and 20°.

Referring to the study of the crack morphology at both tips, the results obtained for the ai = 20° case are similar to those obtained for the ai = 10° case, Section 5.1: contact is detected at the upper crack tip for all stages of growth considered, whereas a ‘bubble’ evolves at the lower crack tip. As was explained in Section 5.2, the detection of the position of the lower crack tip where a physically relevant contact zone appears is a key point for the assessment of the end of the unstable propagation of the interface crack. In this sense a study of the appearance and evolution of the contact zone at the lower crack tip has been carried out for three different initial lengths: ai = 10°, 14° and 20°, the results being presented jointly in Fig. 14. The values of the length of the contact zone shown in Fig. 14 have been estimated interpolating between the results given by the BEM code for the two nodes located at the extremes of the latest element in contact, in the same way as proposed by Blázquez at al. [40], in order to increase the accuracy of the value obtained with the mesh employed. In any case, the size of the last element of the contact zone existing at each position of the crack tip considered is represented by vertical segments included in Fig. 14. The results presented in Fig. 14 show that although the appearance of the contact zone seems to be slightly favoured by a smaller value of ai, there are no relevant differences in the length of the contact zone detected versus the lower crack tip position for the three cases considered. In any case the effect of the elements size undoubtly affects the results obtained. Thus, focusing on the ai = 20° case, it is important to notice that, although the end of the unstable growth was fixed from the energetic analysis at a ffi 200°, the ‘bubble’ at the lower crack tip is still present at this position. In this situation it is plausible to consider a later ending of the interface crack growth up to the position where a contact zone of physically relevant size develops and the effect of friction will unavoidably stop the crack growth, see Section 2.3. 6 initial length: 20º

5

initial length: 14º initial lenght: 10º

ρc (º)

4

3

2

1

0 195

200

205

210

Position of the crack tip, α (º)

215

Fig. 14. Contact zone length at the lower crack tip for the cases ai = 10°, 14° and 20°.

220

4099

E. Correa et al. / Engineering Fracture Mechanics 75 (2008) 4085–4103

Finally, the combined analysis of the results provided by the studies presented in this section makes it possible to conclude that, although the energetic study matches the end of the unstable growth of the crack with a difference of about 6° between the bounds of the range of plausible ai considered, the agreement in the appearance and evolution of the contact zone at the lower crack tip proves a low sensitivity of the propagation of the interface crack to the initial length considered. 5.3.2. Lower crack tip propagation assumption With reference to the initial consideration of the crack enlargement occurring only at the lower tip, two studies will be performed comparing the energetic possibilities of propagation at both crack tips for two different types of growth. The first study is related to the situation studied so far, i.e. a crack propagation occurring only at the lower crack tip, whereas the upper crack tip is fixed at a = 130°. The second study refers to the symmetrical propagation of a crack centred at a = 135°. Concerning the first study, the energetic situation at both crack tips is presented in Fig. 15 for the growth under consideration, i.e. upper crack tip fixed at a = 130° and lower crack tip advancing along the interface. The results are represented in terms of the dimensionless critical load factor, defined as pffiffiffiffiffiffiffiffiffiffiffi rc =r0 ¼ Gc =G; ð17Þ corresponding to both crack tips for each crack size considered. The critical load factor then represents the fraction of the load that is needed to promote crack propagation for each debond considered. In this sense, a value of rc lower than r0 means possibility of growth, whereas a value of rc higher than r0 means no possibility of growth. Referring to Gc, appearing in Eq. (17), it is important to recall that Gc = Gc(wK), it being possible for wK to be different at each crack tip. From the evolutions observed in Fig. 15, it can be deduced that once the growth of the original debond of 10° length is assumed, it is possible for the crack to grow in both directions. In fact, according to this figure, the possibilities of growth are approximately the same for both tips for a crack size lower than 16°, although for debonds larger than 16° growth from the lower crack tip seems more probable. Considering the aforementioned results it appears fundamental to study quantitatively the initiation of growth of a symmetrically situated initial debond with respect to a = 135°. The results of this study are presented in Fig. 16, where G of both crack tips appears versus the size of the interface crack. In this figure it can be observed that the values of G coincide for both crack tips up to a crack length of 45°, where G of the lower crack tip starts to be higher than G of the upper one. In any case, in order to be able to make growth predictions, it is necessary to perform a G comparison between both crack tips in the same way as was previously done for the case of the interface crack growing from its lower tip, Fig. 15. The results for the initial debond situated symmetrically with respect to a = 135° are shown in Fig. 17 where again, once it is assumed that the original debond is able to grow, the possibility of propagation is assured for both crack tips, although now the lower crack tip seems to be more favourable to propagation than the upper one. Recall that although, from an energetic point of view, the possibility of propagation exists for both crack tips, the reasoning for the crack growth prediction requires the consideration of the morphology of the crack at both tips. In this sense, as was explained in Section 5.1, the upper crack tip is completely closed for all situations presented in this paper, whereas a ‘bubble’ exists at the lower crack tip for situations considered in Figs. 15–17. Then, the combination of the conclusions de-

1.2

Upper crack tip Lower crack tip

1

(Gc /G )

0.5

0.8

0.6

0.4

0.2

0 10

15

20

25

30

35

40

45

50

Interface crack size (º) Fig. 15. Comparison between both crack tips of the dimensionless critical load factor associated to the interface crack, considering only growth at the lower crack tip.

4100

E. Correa et al. / Engineering Fracture Mechanics 75 (2008) 4085–4103

0.09

Upper crack tip

0.08

Lower crack tip 0.07

G /G0

0.06 0.05 0.04 0.03 0.02 0.01 0.00 10

15

20

25

30

35

40

45

50

Interface crack size (º) Fig. 16. ERR comparison between the upper and the lower crack tip associated to the interface crack centred at a = 135°.

1.2

Upper crack tip

1

Lower crack tip

(Gc /G )

0.5

0.8

0.6

0.4

0.2

0 10

15

20

25

30

35

40

45

50

Interface crack size (º) Fig. 17. Comparison between both crack tips of the dimensionless critical load factor associated to the interface crack centred at a = 135°.

rived from the information presented in Figs. 15–17 and the arguments relative to the morphology of the interface crack support the lower tip growth considered from the beginning of this section. To sum up the results of this section, it can be concluded that once a small debond (about 10° length centred at the position defined by 135°) is originated by the shear stresses and the critical value of Gc is reached at the lower tip, it will grow unstably as far as a lower crack tip position close to 206°. After this situation, represented in Fig. 18, there are two possibilities: either the crack continues growing along the interface, or the crack kinks and penetrates into the matrix. A further growth of the crack along the interface would necessarily imply an increment of the load. 6. Discussion of the relevance of the present study for dilute and densely packed unidirectional plies Although in the present work the initiation and growth of the fibre–matrix debond under transverse compressive load has been studied in the case of a single fibre embedded by an infinite matrix, it is expected that the conclusions obtained here will also be valid for plies with fibre packing of low density having a sufficient distance between neighbouring fibres, which prevents a strong mutual influence on the interface stresses. Let the relative distance between two neighbouring fibres be

E. Correa et al. / Engineering Fracture Mechanics 75 (2008) 4085–4103

3

σ0

large contact zone

4101

130º

2

initiation of a contact zone

FIBRE 206º

MATRIX

σ0 Fig. 18. Final configuration of the debonded crack.

denoted as d = fibre distance/fibre radius and the inclination angle of the line joining the two fibre centres measured from the external load direction be denoted as h. According to the studies of configurations of multiple fibres embedded by an infinite matrix (e.g., two fibres by Zou and Li [36], four fibres by Crouch and Mogilevskaya [37] and an infinite number of periodically distributed fibres by Isida and Igawa [38]) it seems that the influence of the neighbouring fibres is significant only below a threshold value for d, whose value appears to be about 2 and larger than 1. Recall that in the regular hexagonal packing d = 2 and 1 correspond, respectively, to 23% and 40% fibre volume content. The characterisation of the fibre–matrix debond initiation and growth is much more complex for plies with fibre packing of high density, and in particular for plies, used in engineering practice, with a random distribution of fibres, where the minimum value of d can be very small. Nevertheless, the fundamental results of an analytical study of the interface stresses in the limit case of two touching fibres (d = 0) by Zou and Li [36] can serve as a guide in making predictions about the relevance of the mechanism of interface damage studied in the present work also for these cases. Two particular configurations, corresponding to h = 0° and 90°, were considered in [36]. The influence of a touching fibre on the normal interface stresses is strong but very localized to the neighbourhood of the contact point, where a stress concentration takes place. Considering an external compressive load, a peak of compressions (usually a much stronger one) or tractions appears in the configuration h = 0° or 90°, respectively. The semiangle of the zone of a relevant influence on the normal stress values situated at the contact point is about 25–30°. The distribution of shear stresses in presence of a touching fibre does not change qualitatively in comparison with the case of a single fibre, the maximum difference between these distributions taking place, surprisingly, at an angle of about 50° measured from the contact point. In view of these results, it is expected that the basic qualitative features of the damage mechanisms studied in the present work could characterize the debond initiation and growth at those sides of the fibres (in the approximated angle a range from 130° up to 210°, or any other analogous angular segment) where there is not a very close neighbouring fibre. It might be thought that, according to the above discussed results due to Zou and Li [36], a different mechanism of initiation of the fibre–matrix debond will possibly take place in configurations with h 90°. Note that in these configurations a tensile stress concentration is expected at the closest interface points, with a value of the concentration factor being over 2 for glass fibre and epoxy resin. Nevertheless, it seems, according to the explanations given with reference to Fig. 2, that the tips of such a debond will be closed and a further growth of the debond might be difficult due to the friction dissipation in the contact zones adjacent to these tips. As a conclusion from the above analysis, it is expected that the interface damage mechanism studied in the present work can be relevant, at least from a qualitative point of view, not only for dilute but also for densely packed unidirectional plies. In this sense, the present study establishes a solid base for further more complex studies including several regularly or randomly distributed fibres and covering also subsequent crack kinking towards the matrix. 7. Conclusions This work focuses on the inter-fibre failure of fibrous composites under compression, whose initiation is detected at the fibre–matrix interface in the form of small debonds. The growth of these debonds, becoming interface cracks, has been ana-

4102

E. Correa et al. / Engineering Fracture Mechanics 75 (2008) 4085–4103

lysed by means of a single fibre BEM model. The analysis of the obtained results has been carried out using parameters of interfacial fracture mechanics. The study is initiated considering the fibre completely bonded to the matrix, finding that the most plausible origin of damage is controlled by the shear stress generated between fibre and matrix, the maximum values being placed at the equivalent (due to symmetry) angles ±45° and ±135°, with respect to the direction normal to the compression applied. Once a small debonding crack appears (in particular one of 10° length and centred at 135° has been assumed in the present analysis), its growth is controlled by interfacial fracture mechanics parameters. This growth has been demonstrated to be unstable up to a debond approximately extending from one crack tip at 130° to the other one at 206°. In addition to the relatively large contact zone always existing at the presumably fixed crack tip at 130°, a new small contact zone appears at the crack tip at 206°, a further growth being only possible under Mode II. At this point crack propagation along the interface is only possible if an increment in the applied load is produced (stable growth). The study of the other possibility of growth consisting in a change in the direction of propagation of the crack, i.e. the kinking of the interface crack into the matrix, will be presented in a forthcoming paper [39]. The conclusions obtained have been deduced by means of the stress state at micromechanical level, using the concepts of interfacial fracture mechanics. In this sense, the most relevant concept in the characterisation of the debond growth from a qualitative point of view has been the dependence of the fracture toughness for the debonding crack on the fracture mode mixity, the functional dependences of the fracture toughness assumed in this study being widely accepted in the literature. It should be stated that the numerical solution presented in this work has a clear physical meaning, as the BEM model used is able to detect any physically relevant contact between crack faces, whereas if no such contact is detected these faces can be considered perfectly traction free from the continuum mechanics point of view. A high accuracy of this solution has been confirmed by the self-checking of its coherence at several instances. We can refer, for example, (i) to an excellent agreement achieved when evaluating the phase angle wK of the SIF at the lower crack tip by means of the traction vector at a point ahead of the crack (wK(r)) and by the VCCT integrals (wK(G)), see Fig. 9, and (ii) to a satisfactory agreement between semianalytic estimations of the contact zone length at the lower crack tip (using numerically computed values of wK) with that observed numerically. The conclusions obtained allow an essential stage of the inter-fibre failure of fibrous composites under compression, the damage evolution at the interface, to be characterised. The final situation of the cracks at the interface, predicted by the obtained results, constitutes the key point for the study of the following steps of the mechanism of failure: the growth through the matrix. It is expected that the conclusions obtained will be representative not only for dilute packed unidirectional plies but also for densely packed ones, to the limited extent described in Section 6. This establishes the fundamental basis for a complete knowledge of the failure mechanism, thus helping to build the necessary tools to control its appearance. Acknowledgements The work was supported by the Spanish Ministry of Education and Science (Projects TRA2005-06764 and TRA2006 08077) and Junta de Andalucía (Projects of Excellence TEP 1207 and TEP-02045). The authors thank Dr. C.G. Dávila (NASA Langley Research Center) for his suggestions and comments and Dr. E. Graciani (University of Seville) whose BEM code has been used. References [1] París F. A study of failure criteria of fibrous composite materials. NASA/CR-2001-210661; 2001. [2] París F, del Caño JC, Varna J. The fibre–matrix interface crack. A numerical analysis using boundary elements. Int J Fract 1996;82:11–29. [3] Varna J, París F, Del Caño JC. The effect of crack-face contact on fiber/matrix debonding in transverse tensile loading. Compos Sci Technol 1997;57:523–32. [4] París F, Correa E, Cañas J. Micromechanical view of failure of the matrix in fibrous composite materials. Compos Sci Technol 2003;63:1041–52. [5] París F, Correa E, Manticˇ V. Kinking of transverse interface cracks between fibre and matrix. J Appl Mech 2007;74:703–16. [6] Correa E, Gamstedt KG, París F, Manticˇ V. Effect of the presence of compression in transverse cycling loading on fibre–matrix debonding in unidirectional composites plies. Compos Pt A: Appl Sci Manuf 2007;38:2260–9. [7] Williams ML. The stress around a fault of crack in dissimilar media. Bull Seismol Soc Am 1959;49:199–204. [8] Rice JR. Elastic fracture mechanics concepts for interfacial cracks. J Appl Mech 1988;55:98–103. [9] Comninou M. The interface crack. J Appl Mech 1977;44:631–6. [10] Dundurs J. Discussion of a paper by D.B. Bogy. J Appl Mech 1969:650–2. [11] Manticˇ V, Blázquez A, Correa E, París F. Analysis of interface cracks with contact in composites by 2D BEM. In: Guagliano M, Aliabadi MH, editors. Fracture and damage of composites. WIT Press; 2006. p. 189–248. [12] Graciani E, Manticˇ V, París F. On the estimation of the first interpenetration point in the open model of interface cracks. Int J Fract 2007;143:287–90. [13] Hills DA, Barber JR. Interface cracks. Int J Mech Sci 1993;35:27–37. [14] Irwin GR. Analysis of stresses and strain near the end of a crack transversing a plate. J Appl Mech 1957;24:361–4. [15] Malyshev BM, Salganik RL. The strength of adhesive joints using the theory of cracks. Int J Fract Mech 1965;1:114–28. [16] Manticˇ V, París F. Relation between SIF and ERR based measures of fracture mode mixity in interface cracks. Int J Fract 2004;130:557–69. [17] Comninou M, Schmueser D. The interface crack in a combined tension–compression and shear field. J Appl Mech 1979;46:345–8. [18] Gautesen AK, Dundurs J. The interface crack under combined loading. J Appl Mech 1988;55:580–6. [19] Leguillon D. Singular stress field at the tip of a closed interface crack. In: Ali Mehmeti F, Nicaise S, Von Below J, editors. Lecture notes in pure and applied mathematics. New York: Marcel Dekker; 2001. p. 121–34. [20] Audoly B. Asymptotic study of the interfacial crack with friction. J Mech Phys Sol 2000;48(9):1851–64. [21] Hutchinson JW, Suo Z. Mixed mode cracking in layered materials. Adv Appl Mech 1992;29:63–191. [22] Evans AG, Rulhe M, Dalgleish BJ, Charalambides PG. The fracture energy of bimaterial interfaces. Metall Trans A 1990;21A:2419–29.

E. Correa et al. / Engineering Fracture Mechanics 75 (2008) 4085–4103

4103

[23] Liechti KM, Chai YS. Asymmetric shielding in interfacial fracture under in-plane shear. J Appl Mech 1992;59:295–304. [24] Banks-Sills L, Ashkenazi D. A note on fracture criteria for interface fracture. Int J Fract 2000;103:177–88. [25] Stringfellow RG, Freund LB. The effect of interfacial friction on the buckle-driven spontaneous delamination of a compressed thin film. Int J Solid Struct 1993;30:1379–95. [26] Sun CT, Qian W. A treatment of interfacial cracks in the presence of friction. Int J Fract 1998;94:371–82. [27] París F, Cañas J. Boundary element method. Fundamentals and applications. Oxford: OUP; 1997. [28] Graciani E. Formulation and implementation of BEM for axisymmetric problems with contact. Application to characterization of the fibre–matrix interface in composites. PhD thesis, University of Seville; 2006 [in Spanish]. [29] Graciani E, Manticˇ V, París F, Blázquez A. Weak formulation of axi-symmetric frictionless contact problems with boundary elements. Application to interface cracks. Comput Struct 2005;83:836–55. [30] París F, Del Caño JC, Varna J. BEM analysis of the contact problem in fibres debonded of a matrix. Effect of curing stresses. In: Kassab A, Brebbia CA, Chopra M, editors. Boundary elements XX. CMP; 1998. p. 145–56. [31] Murakami Y. Stress intensity factor handbook. Oxford: Pergamon Press; 1988. [32] Toya M. A crack along the interface of a circular inclusion embedded in an infinite solid. J Mech Phys Solid 1974;22:325–48. [33] Goodier JN. Concentration of stress around spherical and cylindrical inclusions and flaws. 1933;55(7):39–44. [34] Budinsky K. Engineering materials, properties and selection. 2nd ed. Reston PC; 1983. [35] Soden PD, Hinton MJ, Kaddour AS. Lamina properties, lay-up configurations and loading conditions for a range of fibre-reinforced composite laminates. Compos Sci Technol 1998;58:1011–122. [36] Blázquez A, Manticˇ V, París F, McCartney LN. Stress state characterization of delamination cracks in [0/90] symmetric laminates by BEM. Int J Solid Struct 2008;45:1632–62. [37] Zou Z, Li S. Stresses in an infinite medium with two similar circular cylindrical inclusions. Acta Mech 2002;156:93–108. [38] Crouch SL, Mogilevskaya S. Loosening of elastic inclusions. Int J Solid Struct 2006;43:1638–68. [39] Isida M, Igawa H. Analysis of a zig-zag array of circular inclusions in a solid under uniaxial tension. Int J Solid Struct 1991;27:1515–35. [40] Correa E, Manticˇ V, París F. A micromechanical view of inter-fibre failure of composite materials under compression transverse to the fibres. Compos Sci Technol, 2008. doi:10.1016/j.compscitech.2008.02.022.