Numerical investigation into failure of laminated composite T-piece specimens under tensile loading

Numerical investigation into failure of laminated composite T-piece specimens under tensile loading

Composites: Part A 43 (2012) 1017–1027 Contents lists available at SciVerse ScienceDirect Composites: Part A journal homepage: www.elsevier.com/loca...

4MB Sizes 1 Downloads 73 Views

Composites: Part A 43 (2012) 1017–1027

Contents lists available at SciVerse ScienceDirect

Composites: Part A journal homepage: www.elsevier.com/locate/compositesa

Numerical investigation into failure of laminated composite T-piece specimens under tensile loading F. Hélénon ⇑, M.R. Wisnom, S.R. Hallett, R.S. Trask Advanced Composites Centre for Innovation and Science, University of Bristol, Queen’s Building, University Walk, Bristol BS8 1TR, UK

a r t i c l e

i n f o

Article history: Received 9 March 2011 Received in revised form 26 January 2012 Accepted 20 February 2012 Available online 26 February 2012 Keywords: C. Finite element analysis (FEA) A. Polymer–matrix composites (PMCs) B. Stress concentrations Cohesive zone elements

a b s t r a c t This paper presents a rigorous numerical investigation into the structural response of a composite laminate T-piece specimen subjected to a mechanical ‘‘pull-off’’ load case. Initially, a linear elastic stress analysis is conducted, showing very high stresses at the free-edge. In a further analysis, special-purpose interface elements are then inserted where appropriate and used to predict both the crack pattern and the load to failure. It is demonstrated that that using realistic cohesive maximum strength values requires a very fine mesh. Reducing the values to ensure initiation occurs leads to conservative and mesh independent predictions and that a suitable choice leads to good correlation with the experimental results. This study also shows that the T-piece failure is controlled by crack propagation. Ó 2012 Elsevier Ltd. All rights reserved.

1. Introduction The present paper deals with the finite element (FE) analysis of a T-cross section structure made of carbon fibre reinforced plastic which is typically encountered in vanes, stiffeners and in other beam structures. Because of their relative complexity, the prediction of both their damage resistant capability and of their failure mechanism is fundamental for a conservative design. The type of T-shaped composite structures addressed in this paper is made of three sections laid up from pre-preg material and a fill-in region at their junction; see Fig. 1. Two1 of them (in blue) are formed at 90° such that the upper vertical parts bonded together form the ‘web’. The radiused corners are referred to as ‘fillets’. The horizontal parts overlap the third section (in green) namely the ‘platform’ or the ‘skin’ to form the ‘flange’. At the junction there is a fill-in region (in red) named ‘deltoid’ or ‘noodle’. This reinforcement block may be filled in with resin, tows or with adhesive to help preserve the joint structural integrity. As discussed in Panigrahi and Pradhan [1], four different types of deltoid profiles can be made during the design process, namely triangular, circular, elliptic or modified elliptic. One of the early reported works on such a T-piece geometry was by Gillespie and Pipes [2] who modelled joints under pull-off loading in a composite spar concept. Joints with or without a horizontal titanium or graphite-epoxy insert within the wingskin were ⇑ Corresponding author. Tel.: +44 117 331 5506; fax: +44 117 954 5666. E-mail address: [email protected] (F. Hélénon). For interpretation of color in Figs. 1–13, the reader is referred to the web version of this article. 1

1359-835X/$ - see front matter Ó 2012 Elsevier Ltd. All rights reserved. doi:10.1016/j.compositesa.2012.02.010

considered as a structural reinforcement. Correlation between experimental load–deflection response and that from a plane stress linear elastic FE analysis showed an excellent agreement. Cope and Pipes [3] designed various spar–wingskin joint concepts with triangular and curved fillets in order to determine the most conservative deltoid profile. Plane stress FE models were used in their investigation and the Tsai-Wu and maximum stress failure criteria were used to identify the delamination onset location. This choice provided bounds to most experimental data, and the authors demonstrated that design concepts with curved joints maximises the overall strength compared to triangular ones. Furthermore, verifications have shown that increasing the radius in the corner increases the overall joint strength. Later, Rispler et al. [4] evaluated the influence of different insert types placed within the resin rich deltoid area by correlating experimental results and plane strain linear elastic FE analyses. They showed that their FE models were able to accurately predict failure initiation location on most specimens regardless of insert type and web thickness. They also concluded that, among several failure criteria, the interlaminar peel failure criterion provided the most accurate and conservative failure initiation load prediction. In an effort to provide better insight to T-shaped structure designers, Kumari and Sinha [5] performed a numerical parametric study using FE analyses of wing T-joints with no filler insert. The joints were modelled as being under pull-out loading and a thick composite shell finite element was developed and employed. By considering two different stacking sequences, they showed that the failure mode is configuration dependent. In addition, they noted that increasing the radius between two fixed points of arc length at the web/skin interface increased the pull-off strength

1018

F. Hélénon et al. / Composites: Part A 43 (2012) 1017–1027

Fig. 1. CFRP T-joint under consideration in this paper. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

because stress concentrations were reduced at this location. This was only valid up to a certain point since a further increase became ineffective. In contrast to the work of Cope and Pipes [3], the authors also claimed that when decreasing the filler area by varying the end-point positions of arc length, the joint strength increases. They state that this was valid only down to a certain area because further decreases would cause fibre kinking during manufacture. To gain better understanding of the joint’s behaviour during service life and knowing that T-joints are prone to delamination due to the lack of through-thickness reinforcement, Rao et al. [6] investigated both experimentally and numerically the effect of moisture on transversely stitched and unstitched composite aircraft wing T-joints subjected to a pull-off force. The experiments were conducted under dry and hygrothermal conditions with a 0.7% saturated moisture concentration. Joints having respectively 2.4 mm and 3.6 mm corner radii were studied. Both linear and non-linear FE analyses undertaken and compared to experimental measurements. The authors found that moisture leads to about 20% reduction in the strength regardless of the joint configuration. Only the joint having a 3.6 mm skin thickness behaved almost linearly up to failure. The authors showed that increasing the skin thickness enhances both the failure initiation and final failure load. They also concluded that transverse stitching enhances the failure load, the failure mode and the failure pattern. More recently, Zimmermann et al. [7] studied ultra-thick composite joints whose thickness exceeds 60 mm, as part of the development of the side stay fitting in a main landing gear. Pieces manufactured entirely with non-crimped fabrics were tested in tension, axial bending and in compression. 3D linear elastic FE approaches using 8-noded stacked composite elements were chosen due to the presence of transverse shear and through-thickness stresses. For tension, high localised peeling stress concentrations were noticed at the free-edge of the joint and were responsible for failure at the upper deltoid tip through the web. For bending, the authors showed that additional significant shear stresses in the fillet contribute to failure. Finally, when loading in compression, the authors have shown that using a stiffer gusset filler made of carbon instead of foam is better since it diverts the load and significantly reduces the shear stresses and the downward displacement of the joint. Predicting failure initiation is important for design, and the prediction of delamination induced damage is also essential for assessing the performance of the joint. This was carried out by by Panigrahi and Pradhan [1] who used a 3D non-linear FE analysis

to predict damage onset and growth in a composite spar–wingskin joint having a modified elliptic deltoid profile fully filled with adhesive. The Tsai-Wu stress based failure criterion was used for predicting damage initiation and VCCT was used to model delamination propagation along a predefined crack path at the spar/wingskin interfacial surface. The authors concluded that the complete failure of the joint was due to delamination propagation initiating from the toe-ends of the interfacial surface. This motivated them to extend this work and to study the influence of different deltoid profiles on the strength of the composite spar wingskin joint [8]. Triangular, circular, elliptical and modified elliptical load couplers were considered. By analysing the stress state and by using the Tsai-Wu stress based criterion to predict delamination onset damage, they made several recommendations for the design of composite spar wingskin joints. Indeed, they have found that circular or elliptical fillet profiles are preferred over triangular ones, which should be avoided. They also found that modified elliptical profiles are the most suitable and an optimum design was identified as being with a ratio of 1.5 between the base width and the height of the deltoid area. The same year, Chen et al. [9] used special-purpose cohesive elements in Abaqus/Standard in order to study the progressive delamination within a braided composite T-piece specimen modelled in plane strain. Elements were inserted only at the critical layer subjected to delamination failure as noticed from experimental observations. Modelling predictions showed that the damage behaviour was in good agreement with experimental observations. In addition, by varying independently the fracture toughness properties and the strength values by ±50%, the authors showed that failure was dominated by mode II dominated mixedmode crack propagation rather than its initiation. In-line with Chen et al’s work, cohesive interface elements are used in the work presented in this paper to conduct a rigorous numerical investigation into the structural response of a laminated composite T-piece specimen under ‘‘pull-off’’ loading. Initially, both 2D and 3D linear elastic FE analyses are carried out and compared. The importance of the free-edge effect, which was missing in most of the reported work summarised above, is also discussed. The transverse maximum principal stress criterion is shown to be relevant in both the 2D and 3D FE approaches, and is used to insert special-purpose interface elements into non-linear FE models, independently from experimental observations to predict both the crack pattern and the failure load. The results are found to be in good agreement with experimental results and it is demonstrated that the T-piece failure mechanism is dominated by crack propagation.

2. Experimental set up The experimental details of the T-piece specimen under consideration in this paper are presented in Trask et al. [10]. A brief summary of the important and relevant results are recounted here for completeness. Specimens were manufactured from Hexcel’s IM7/ 8552 unidirectional carbon/epoxy pre-preg tape in three separately laid-up sections: two curved and one straight with a fill-in deltoid region made from 90° unidirectional tows. Thermo-elastic properties are given in Table 1. Each specimen is 20 mm wide and has a 5 mm external radius corner. An overview of the stacking sequence around the lower right deltoid tip is provided in Fig. 2. The nominal ply thickness is 0.127 mm. The flange (horizontal base) is 5.334 mm thick with a (60/0/60/0)3S layup for the platform plus a (±452/07/904/03) lay-up which has come down from the web. At both of its ends, the flange is in contact on its upper surface with two rollers having an 80 mm span. They prevent the piece from moving vertically when pulled upwards by the web. The latter is 4.572 mm thick with a (±452/07/904/03)S lay-up. The

F. Hélénon et al. / Composites: Part A 43 (2012) 1017–1027

web is firmly clamped into jaws at 40 mm from the underside of the flange, and subjected to a controlled vertical displacement. The T-piece reaction force and cross-head displacement is then measured and the load–displacement curve plotted, see Fig. 3. A linear behaviour is noticed up to failure, when a sudden load drop takes place. For specimens having a 24 mm2 deltoid area, the average failure load from five tests has been reported as 1330 N with a 9 % coefficient of variation [10]. Fig. 4 highlights what generally occurs at failure when the Tpiece is under pull-off loading. It shows a crack pattern to have initiated within the 904 block and propagated at both the 07/904 and 904/03 interfaces. This localisation of damage indicates that the 904 block is the most sensitive part of the specimen. As reported in [10], other similar T-piece specimens exhibited cracks occurring mainly within the deltoid region. However, microscopic analyses have suggested that such a pattern is due to the quality variation in manufacture. Indeed, a higher number of voids were present in the very thick block of deltoid material, which has a direct influence on damage initiation [4,11]. 3. Damage onset and predefined crack pattern prediction In order to predict the damage onset locations, a linear elastic approach is used. Several researchers have already used such an approach successfully for similar purposes [2–5,12]. In the rest of this section, results from plane stress and 3D approaches are therefore discussed. 3.1. Plane stress approach A plane stress approach is firstly undertaken, motivated by the good predictions made by Gillespie and Pipes [2] and Cope and Pipes [3] among others. Such a modelling assumption is reasonable since the specimen is 20 mm wide, provided free-edge effects are not significant. A load of 33.25 N/mm width is prescribed at the top of the web, consistent with the 1330 N average experimental failure load. Since the T-piece geometry and loading is symmetric, only half of the geometry, as given in Fig. 2d, is meshed via Abaqus/ CAE; see Fig. 5. Symmetric boundary conditions (vertical row of black triangles) are used along the middle plane of the web; see Fig. 2b. The last node at the top end of the flange is constrained vertically to account for the presence of the roller during the pull-off loading; the remaining flange portion outside the 80 mm span is neglected in the FE model because it does not contribute to the state of stresses at the flange/web junction. As it can be noticed in Figs. 2d and 4, the vertices of the deltoid region have been redrawn finishing slightly upstream of the actual geometry to avoid very distorted elements in those regions. The small resin-rich pockets at the deltoid tips have not been modelled. Although they may affect the stress state locally, this has a negligible effect on the stress state within the surrounding material. Around the deltoid, 2 elements per ply are used to model the (±452/07/904/03)S curved region and 1 element per ply is used to model the (60/0/60/0)3S platform. This gives a mesh with 12,179 bi-quadratic elements of type CPS8 everywhere except at the deltoid tips where two triangular quadratic elements of type CPS6 are introduced. Equivalent orthotropic thermo-elastic material properties are used to model the (60/0/60/0)3S platform and the ±452 block since out-of-plane shear stresses and twisting cannot take place in 2D models. Table 1 Thermo–elastic properties for IM7/8552 [18]. Subscripts ‘l’ and ‘t’ denote longitudinal and transverse directions respectively. El = 161 GPa Et = 11.38 GPa

mlt = 0.32 mtt = 0.436

Gtt = 3.98 GPa Glt = 5.17 GPa

all = 0 /°C1 att = 30 /°C1

1019

The calculations are broken down into two distinct steps. The first models the cure phase with a 160 °C variation of temperature applied to the whole specimen. During this step, only the symmetric boundary conditions are applied. One node is however fully constrained to prevent the model from moving rigidly. The second step focuses on mechanical stresses generated when applying the 33.25 N/mm width failure load in addition to the residual stresses. From this decomposition, it is found that the curved deformation of the flange is predicted but none of the residual stress tensor components are high enough to cause premature failure onset damage. Focussing on the second step which accounts for both the superimposed residual and mechanical stresses, it is found that the most highly loaded regions are the 904 block and the deltoid regions; see Fig. 6. The stress field is fairly difficult to interpret because the normal (S11), the transverse (S22) and shear stresses (S12) are all present together. They have therefore been combined into the maximum principal stress perpendicular to the fibre direction, which captures the contribution of the different components affecting matrix dominated failure. Such an approach has been found to give satisfactory results by researchers such as Minguet and O’Brien [13] when studying composite panels reinforced with stringers. Fig. 6 shows the most sensitive parts determined from the plot of this stress tensor component i.e. the deltoid and the curved region of the 904 block. The highest stress value (93 MPa) is found in the curved region of the 904 block. By assuming the transverse tensile strength to be 92 MPa [14], failure would be expected to occur at that location for the applied load of 33.25 N/mm width. This conclusion is in agreement with the experimental observations in Fig. 4 showing the presence of a transverse matrix crack at the same location. Moreover, the plot of the principal vector associated with the maximum principal stress component at each integration point shows the direction that the crack should follow once initiated from the most highly loaded point. A good correlation is found again by comparing the predicted direction with the angle of the crack within the 904 block. Once fully formed through the 904 block, the presence of fibres means that the crack has no other choice than propagating at the 07/904 and 904/03 bimaterial interfaces, as observed experimentally.

3.2. 3D approach In the previous section, the plane stress approach has enabled the damage onset location to be determined in agreement with experimental observations with the assumed transverse tensile strength of 92 MPa. This value is however lower than other reported values for IM7/8552 well above 100 MPa [15–17]. Furthermore, other cracks reported within the deltoid region [10] were not predicted and one of the main sources of delamination, namely the free-edge effect, was missing. In order to account for those aspects, a 3D model is developed. An identical block homogenisation approach to the one used in the 2D model is used. This choice is appropriate since no free-edge effect is expected for either the (60/0/60/0)3S platform or the ±452 block. This enables a quarter-geometry to be used as shown in Fig. 7, which reduces considerably the computational costs. Symmetric boundary conditions are used on the Y–Z and X–Y planes. The support (roller in Fig. 2b) is taken into account by preventing all the nodes at the top-end of the flange from moving vertically. The top of the web is loaded with a 14.545 MPa uniform normal stress distribution consistent with the 1330 N average failure load [10]. 70,246 hexahedral and 44 triangular prism elements of type C3D20 and C3D15 respectively are used. Globally, one element per ply is used except for the platform region where there is one element per 2 plies. Across the width, the free-edge is more refined in order to guarantee a good representation of the high stresses. From the free-edge

1020

F. Hélénon et al. / Composites: Part A 43 (2012) 1017–1027

(a)

(b)

(c)

(d)

Fig. 2. (a) T-piece experimental testing configuration, (b) schematic illustration of the test setup, (c) microscopy section detailing fibre stacking sequence around the lower right deltoid tip fill-in with resin, and (d) schematic illustration of the nominal T-piece geometry showing deltoid tips finishing slightly upstream of the actual position (red circles) to avoid extremely distorted elements when meshing. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 3. Experimental load–crosshead displacements curves. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

up to the plane of symmetry, 10 elements 0.2 mm wide then 12 elements starting from 0.23 mm up to 1.4 mm are used. This mesh will be denoted as being ‘‘Mesh B’’ in the rest of this paper. As in the plane stress approach, the calculations are split into two steps to analyse each loading condition’s contribution to the level of generated stresses. It is therefore found that, even with the free-edge contribution, none of the residual stresses are high enough to cause premature failure. Therefore, only the subsequent case with combined residual stress and pull-off loading phase is discussed hereafter. This is performed by plotting the most rele-

Fig. 4. Typical crack pattern (highlighted in red) occurring in some of the experimental tests (recorded failure load here: 1.26 kN). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

vant stress tensor components within the most sensitive blocks of material identified; see Fig. 8. It can be noticed that the highest values have been computed at the free-edge of the 904 block in the corner of the fillet radius. Drawing direct conclusions regarding the failure mode from those stress fields is challenging. Combining them into some form of stress-based criterion such as Tsai-Hill or LaRC would enable the

1021

F. Hélénon et al. / Composites: Part A 43 (2012) 1017–1027

radius corner are observed. Moreover, the well-known free-edge effect has induced a very high localised stress concentration at the free-edge of the 904 block. Such 3D FE stresses imply that the crack observed experimentally initiates from the free-edge of the specimen. The predicted failure load is also expected to be much lower than the measured one, even if a higher material transverse tensile strength such as 111 MPa [16] is assumed. 4. Delamination damage growth prediction using interface elements

Fig. 5. Schematic illustration of the plane stress approach for the T-piece specimen and overview of the 2D mesh in the deltoid zone. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

onset damage location to be identified but would not provide information regarding the initial crack propagation direction and, therefore, the expected crack pattern. This is why the maximum principal stress perpendicular to the fibre is chosen again even for this 3D FE analysis. By plotting in Abaqus/Viewer the vectors associated with the maximum principal stresses and by looking at their orientation from the T-piece lateral side (X–Y plane), it is possible to predict the orientation of initial crack propagation, see Fig. 9 which can be compared to the prediction made in Fig. 6d. Similar stress concentration locations within the 904 block

Special-purpose interface elements developed at the University of Bristol [18] are now considered in order to extend the previous FE analyses to both the crack pattern formation and average failure load predictions. They take the form of solid hexahedral elements with a small initial thickness (0.01 mm) and are governed by the single three-dimensional map in Fig. 10 representing the normal opening mode (mode I) on the 0–rII–dI plane, and the transverse shear mode (mode II) on the 0–rII–dII plane. The triangles 0– f rMax –dfI and 0–rMax I II –dII are the bi-linear responses in the pure opening and pure shear modes respectively. Any point on the 0– dI–dII plane represents a mixed-mode relative displacement. The mixed-mode damage onset displacement, dem , and interfacial strength, rMax m , are calculated using the quadratic damage onset criterion.

sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  2  ffi maxðrI ; 0Þ rII 2 ¼1 þ Max Max

rI

rII

ð1Þ

The mixed-mode failure displacement corresponding to complete decohesion is calculated using the following power-law failure criterion:

GI GII þ ¼1 GIc GIIc

ð2Þ

where, GIc and GIIc are the critical energy release rates for pure mode I (opening) and pure mode II (shear), respectively. Eq. (2) allows the fully debonded locus, represented by the relative displacement corresponding to complete interface failure, dfm , to be determined.

(a)

(b)

(c)

(d)

Fig. 6. Combined residual and mechanical FE stress fields within deltoid area and curved region of 904 block for a 32.25 N/mm width applied pull-off load. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

1022

F. Hélénon et al. / Composites: Part A 43 (2012) 1017–1027

Fig. 7. Overview of the quarter FE model used for preliminary analysis (Mesh B). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 8. Relevant stress field components within the critical regions (deltoid and 904 block) at the experimental failure load. The black arrows indicate the considered directions for the stress components. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

In order to use the interface elements presented above, the same Abaqus FE model given in Fig. 7 is considered but now converted into an LS-Dyna FE model having 8-noded linear elements instead of 20-noded ones. The mesh is thus effectively coarser, and the issue of mesh refinement is considered later. The calculations are carried out by modelling the cure phase first, by applying a 160 °C reduction of temperature to the whole FE model. Afterwards, a vertical displacement at the top-end of the web is prescribed and the nodes at the upper-end surface of the flange are constrained vertically during the pull-off phase. This induces a reaction force from the specimen which can be measured either at the top-end of the web or at the end of the flange.

4.1. Importance of the predefined transverse crack path Two different crack patterns meshed with interface elements are tested in Mesh B, see Fig. 11. In all cases, a 0.01 mm interface element thickness is introduced, which is the typical resin-rich interface thickness value between plies [19]. Such a value is sufficiently small not to affect the surrounded stresses within the laminates. The first case accounts for interface elements inserted at the 07/ 904, the 904/03 and the 03/deltoid interfaces whose element lengths in the hoop direction are about 0.25 mm. This leads to a mesh, denoted as Mesh B1, having 61,974 linear elements in total;

F. Hélénon et al. / Composites: Part A 43 (2012) 1017–1027

Fig. 9. Transverse maximum principal stress field within the deltoid and 904 block and prediction of the location and initial direction of the crack propagation. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 10. Illustration of the mixed-mode response in interface elements. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

see Fig. 11a. The second case accounts for interface elements inserted only at the 07/904 and the 904/03 interfaces as well as for a transverse pre-embedded potential matrix crack through the 904 block placed perpendicularly at the location of peak transverse maximum principal stress, from the linear elastic FE analyses. At the 07/904 and the 904/03 interfaces, the elements lengths in the hoop direction are about 0.1 mm close to the transverse predefined matrix crack tip up to about 0.29 mm close to both ends of the fillet. Those for the transverse predefined matrix crack are about 0.13 mm long. This mesh is denoted as Mesh B2 having 78,298 linear elements in total. Fig. 11b illustrates where the interface elements are inserted in this second case. The IM7/8552 resin-dominated properties used to model the interface element behaviour are given in Table 2 in which the strength values, S and S, and fracture toughness values, GIc and GIIc, have been extracted from [18,20], and the interface stiffnesses, kI and kII, have been selected following [21]. Even though there is evidence that fracture toughness values may be different for a 0°/90° and a 90°/90° interface [22], similar values are used for all the interface elements regardless of their insertion location. Such an assumption is common in practice and enables good predictions to be made [19].

1023

By post-processing the results for both Mesh B1 and Mesh B2, it is found that the T-piece behaves linearly during the pull-off phase up to failure even with interface elements inserted, see Fig. 12. An overall linear response is observed in the traction-displacement curve before occurrence of a sudden load drop when the specimen fails completely. This is in agreement with experimental observations reported by Trask et al. [10] and shown in Fig. 3 and is consistent with behaviours reported by other researchers on other carbon/epoxy laminated T-joints [2–5,9,12]. The predicted stiffness is however slightly higher than the experimental one, but this is typical for numerical models ignoring the compliance of the test fixture when compared against measured cross-head displacements. Regarding Mesh B1, failure initiation and immediate unstable delamination failure occurs only at the 07/904 interface once the first interface element fails. The crack propagates simultaneously in the hoop direction and across the width. This happens when the reaction force reaches 2991 N, a value far above the expected 1330 N average experimental failure load. From Mesh B2 accounting for the transverse predefined matrix crack in the fillet, failure occurs at this location first, then propagates almost instantaneously at both the 07/904 and 904/03 interfaces. Overall, unstable crack growth is observed immediately after initiation, once the first interface element fails, when the load reaches 1655 N, which is still above the 1330 N expected failure load. As in Mesh B1, the crack grows simultaneously across the width and in the hoop direction, giving a crack pattern similar to the one observed experimentally in Fig. 4. Such results show the importance of the transverse matrix crack in the failure mechanism. Indeed, in Mesh B1 with no potential transverse matrix crack, the through-thickness and the interlaminar shear stresses are not high enough to fulfill the quadratic failure onset criterion, see Eq. (1), within the interface elements when the prescribed displacement leading to failure in Mesh B2 is attained. This can be understood by looking at the stress fields plotted in Fig. 8. Even though edge stresses are higher there than those computed in LS-Dyna (because a quadratic FE mesh was used) combining the through-thickness stress S22 with the interlaminar shear stresses S12 and S23 (by using the Abaqus/Viewer tool to calculate a failure index based on the quadratic stress interaction criterion, see Eq. (1)) gives a maximum value of 0.79 at the free-edge, which is far below the value of 1 required to cause onset of failure [18]. Furthermore having instantaneous unstable crack propagation in both cases is due to a large amount of energy of deformation stored in the specimen before its complete failure. Following initiation, high energy release rates fulfill straight away the linear energetic failure criterion for crack propagation, see Eq. (2), once failure initiation occurs. However, when Mesh B2 is used, the most loaded region i.e. the predefined matrix crack path is under mode I opening because it is loaded by the maximum principal stresses. It can therefore be expected to get a lower failure load prediction with this mesh. In this case, the maximum failure index at the free-edge, based on the transverse maximum principal stress component, is equal to 1.2 in the Abaqus FE model. 4.2. Influence of the free-edge mesh refinement Since the obtained FE edge stresses are mesh-dependent, the influence of local mesh refinement is now assessed by considering together Mesh A, Mesh B and Mesh C; see Fig. 13. In each of these meshes, a predefined crack pattern accounting for the transverse matrix crack within the 904 block, i.e. the one used in Mesh B2, is accounted for. All of them have the same lateral mesh definition extruded across the width; see Fig. 11b. However, Mesh A is coarser than Mesh B across the width, built with element sizes starting from 0.5 mm at the free-edge up to 1.38 mm at the X–Y plane of symmetry. Mesh C, finer at the free-edge than all the others, is sim-

1024

F. Hélénon et al. / Composites: Part A 43 (2012) 1017–1027

(a)

(b)

Fig. 11. Overview of the FE meshes from the lateral side with predefined crack paths: (a) Mesh B1 with interlaminar crack paths only (in red). (b) Mesh B2 with interlaminar predefined crack path (in red) and with transverse crack within 904 block (in red). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Table 2 Resin dominated properties for IM7/8552 [18,20,21].

rS (MPa) 111

sS (MPa) 130

GIc (N/mm) 0.2

GIIc (N/mm) 1

kI (N/mm3) 6

10

kII (N/mm3) 106

ilar to Mesh B apart from the first row of elements at the free-edge subdivided into three rows 0.067 mm wide. Irrespective of the FE model, it is found that unstable failure propagation occurs once the first interface element fails, immediately after initiation at the free-edge, but at different loads; see right column of Table 3, ignoring the results for other values of strength parameters for now. Mesh C (refined) shows the earliest appearance of damage, followed by Mesh B2 and then Mesh A (coarse). This is understandable since higher edge stresses are generated when increasing the edge mesh density [23]. This enables the stress-based failure initiation criterion [18] to be fulfilled earlier in the refined mesh than in the coarse mesh. A direct consequence is that the finer the mesh, the lower the predicted failure load, although the predictions are all higher than the expected 1330 N average experimental value. This is because all the meshes

are in fact too coarse, and so failure is controlled by initiation rather than propagation. If the edge refinement was continued, a further decrease of the predicted failure load would be expected, which might approach the 1330 N experimental value. There would, however, be a significant increase in computation time. More elements would be used in the FE model and a smaller critical time step would be required due to the smaller element size at the free-edge. It is now of significance importance to discuss the numerical cohesive zone length and the relevance of the chosen mesh refinement and strength parameters used so far in this work. Indeed, as discussed in [19,24,25], a mesh has to be fine enough to have several elements in the cohesive zone to give accurate analysis during the crack propagation regime. To help choose the suitable mesh refinement, several formulas estimating the cohesive zone length for mode I and mode II loading cases have been proposed in the literature assuming plane stress conditions. They all have the following form: I

=

lcz ¼ MEI

GIc ðrMax Þ2 I

II

=

and lcz ¼ MEII

GIIc 2 ðrMax II Þ

ð3Þ

where M is a scaling factor, GIc, GIIc are respectively the mode I and mode II fracture toughness values, and rMax , rMax are the maximum I II

Fig. 12. Typical load–displacement curves obtained from the FE models including interface elements with final crack pattern shown. The Mesh B2 curve (with failure load of 1655 N) was offset by +0.01 mm for the sake of clarity. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

F. Hélénon et al. / Composites: Part A 43 (2012) 1017–1027

1025

Fig. 13. Overview across the width of the three mesh refinements defined for the quarter T-piece geometry. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

4.3. Influence of the resin dominated properties Table 3 Predicted failure loads against average experimental value as a function of the FE model and of the selected pairs of maximum stress values considered. Mesh A, Mesh B2, Mesh C and the slice model have the same lateral mesh definition extruded across the width (with potential transverse matrix crack). No transverse crack is introduced in Mesh B1. Pairs of maximum stress values ðrMax ,rMax I II )

FE models

Quarter model

Mesh Mesh Mesh Mesh

Slice model Experimental tests (average value)

A B1 B2 C

(30,60)

(60,90)

(111,130)

966 – 965 965 1009 1330

1274 – 1262 1260 1358

1840 2991 1655 1586 2113

cohesive stresses, equal to the interfacial strengths already given in = = Table 2, and where EI , EII are elastic constants defined in [25]. There is considerable debate regarding the choice of M and the empirical value of 0.5 is used here, as recommended by Harper and Hallett I II [25]. This gives lcz ¼ 0.11 mm and lcz ¼ 1.55 mm and, as a consequence, the interface element length should be ideally below half of the lowest value i.e. 0.05 mm. This condition is not met in any of the meshes. For Mesh B1 the interface element sizes are all over 0.2 mm and for Mesh A, Mesh B2 and Mesh C, the smallest ones are about 0.12 mm within the transverse crack. Therefore, the T-piece FE mesh would need to be significantly refined through the width in order to fulfil this cohesive zone length condition, which is not practical.

In order to overcome the issue of the need of an extreme mesh refinement, a study was undertaken to reduce the strength values given originally in Table 2. As suggested by Turon et al. [19], this choice will cause the numerical cohesive zone length to increase and should allow early failure initiation, without having quasiinstantaneous unstable crack propagation. Therefore, two new pairs of tensile transverse and shear strength values, (rMax ,rMax I II ), are now tested, namely (60,90) and (30,60). Such values were already used in other contexts such as in [26–28] to model simple coupon specimens. Here, they are re-used in this work to evaluate their effect on a more complex model. From Eq. (1) the recommended numerical cohesive element lengths are 0.39 mm and 1.54 mm for each of the pairs of strength values respectively. As a consequence, simulations made from the chosen meshes should be a priori sufficiently accurate during the crack propagation regime, except for Mesh A with the (60,90) strengths since the minimum element length is 0.5 mm. After having run the three FE models presented in Fig. 13 using the (60,90) pair, it is found that damage occurs at the free-edge of the transverse predefined matrix crack path at approximately 60% of the final load. It then growths until the most loaded interface element fails. The crack created there then propagates stably until complete failure of the transverse predefined crack (99.5% final load). Slight delamination crack propagation appears then at the free-edge of the predefined crack paths introduced at both the 07/904 and the 904/03 interfaces. Finally, unstable crack propagation occurs from those delamination locations resulting in complete failure of the specimen. This causes a load drop on the

1026

F. Hélénon et al. / Composites: Part A 43 (2012) 1017–1027

force–displacement curve which remained linear until this event. The highest recorded load value at failure is given in Table 3 according to the FE model considered. It can be seen that all three levels of mesh refinement considered now give similar values, which are within about 5% of the experimental results. When the strength values are further reduced down to the (30,60) pair, damage appears during the cure regardless of the mesh used (Fig. 13). This induces an early stable crack propagation occurring just after having begun to load the specimen without affecting the predicted stiffness which remained similar to the ones obtained with the previous (111,130) and (60,90) pairs. As a result, an early complete failure of the specimen is observed and the corresponding load is reported in Table 3, again shown to be independent of the mesh size used. This mechanism was not observed in the tests and shows that reducing the strength parameters too far can cause premature failure or a change in mechanism. It can be concluded that the Tpiece specimen fails due to propagation of delamination initiated at the 904 block. Panigrahi and Pradhan [1] and Chen et al. [9] have found similar conclusions regarding failure by crack propagation when studying other T-piece specimens loaded up to failure. 4.4. Comparison with a slice model Since independence of mesh size across the width on the predicted failure load using reduced strength values has been obtained, it is interesting to investigate whether or not a 2D approximation, which does not take account of the free-edge, would give a similar failure load prediction. In this regard, a model with one element across the width is considered, i.e. a 0.2 mm slice model built with the same lateral mesh definition as the one used previously in Mesh A, Mesh B2 and Mesh C; see Fig. 11b. Such a thickness is chosen so that the aspect ratio of any element is acceptable and because 2D elements, which would have been more appropriate for this approach, are not available in LS-Dyna. The boundary conditions are all similar to those previously used, see Fig. 7, apart from the free-edge which is constrained by using MPCs to couple all the displacements in the Z-direction to achieve a generalised plane strain boundary condition. The results obtained from this model are summarised in Table 3. It can be noticed that using the (111,130) pair of cohesive strength values leads to an extremely high failure load (2113 N) because initiation cannot occur at the expected failure load of 1330 N with this level of mesh refinement. However when using (30,60), a too conservative prediction is found (1009 N) since the predicted failure mode changes. As with the quarter T-piece FE models, these strength values led to an early onset of damage appearing during the cure phase and extensive stable delamination early on during the loading. This premature failure no longer occurs when using the (60,90) pair, and the mesh is now sufficiently fine to allow initiation. Damage occurs at the transverse predefined matrix crack during the mechanical loading phase for a prescribed displacement similar to the one observed using the 3D quarter Tpiece FE models. Once the interface elements in the transverse crack path are fully damaged, unstable delamination failure then occurs once the energy release rate is sufficient to fulfil the linear damage propagation criterion. A very close agreement (1358 N) is therefore found with the experimental failure load, with the value slightly higher than those obtained from 3D models. 5. Conclusion In this paper, a rigorous numerical investigation into the structural response of a laminated composite T-piece specimen under pull-off loading was carried out. This study started with predictions of the expected crack pattern using the transverse maximum prin-

cipal stress criterion, which was in good agreement with experimental observations. Special-purpose interface elements were then inserted where necessary within the T-piece FE mesh in order to capture the full damage development and predict specimen failure. It was demonstrated that using realistic strength values in the model requires an extremely fine mesh to allow initiation to occur. Using reduced strength values, good failure load predictions were obtained, independent of the mesh refinement at the free-edge of the specimen. However reducing the values too far caused a change in behaviour, with premature failure during the cure. Finally, it was shown that the T-piece specimen fails by delamination propagation and not at the onset of damage since the free edge stress value varies with mesh refinement level in the (60,90) case yet the predicted failure load remains constant. Suitable insertion of through-thickness reinforcement such as stitching or Z-pinning could arrest or even eliminate delamination induced damage and therefore increase the T-piece loading capability. Acknowledgement The authors would like to acknowledge the support of RollsRoyce Plc in funding this work through the composites University Technology Centre (UTC) at the University of Bristol, UK. References [1] Panigrahi SK, Pradhan B. Delamination damage analyses of FRP composite spar wingskin joints with modified elliptical adhesive load coupler profile. Appl Compos Mater 2008;15(4–6):189–205. [2] Gillespie JW, Pipes RB. Behaviour of integral composite joints – finite element and experimental evaluation. J Compos Mater 1978;12:408–21. [3] Cope RD, Pipes RB. Design of the composite spar–wingskin joint. Composites 1982;13(1):47–53. [4] Rispler AR, Steven GP, Tong L. Failure analysis of composite T-joints including inserts. J Reinf Plast Compos 1997;16(18):1642–58. [5] Kumari S, Sinha PK. Finite element analysis of wing T-joints. J Reinf Plast Compos 2002;21(17):1561–85. [6] Rao VVS, Krishna Veni K, Sinha PK. Behaviour of composite wing T-joints in hygrothermal environments. Aircraft Eng Aerospace Technol 2004;76(4):404–13. [7] Zimmermann K, Zenkert D, Siemetzki M. Testing and analysis of ultra thick composites. Composites Part B 2010;41:326–36. [8] Panigrahi SK, Pradhan B. Development of load coupler profiles of spar wingskin joints with improved performance for integral structural construction of aircraft wings. J Reinf Plast Compos 2009;28(6):657–73. [9] Chen J, Ravey E, Hallett S, Wisnom M, Grassi M. Prediction of delamination in braided composite T-piece specimens. Compos Sci Technol 2009;69(14): 2363–7. [10] Trask RS, Hallett SR, Helenon F, Wisnom MR. Influence of process induced defects on the failure of composite T-joints specimens. Composites Part A 2012;43(4):748–57. [11] Wisnom MR. Statistical aspects of failure of fibre-reinforced composites. Proc Inst Mech Eng. Part G: J Aerospace Eng 1998;212(3):189–92. [12] Hill GFJ, Wisnom MR, Jones MI. Failure prediction of composite T-piece specimens. Proceedings of the 5th international conference on deformation and fracture of composites; 1999. p. 10. [13] Minguet PJ, O’Brien TK. Analysis of test methods for characterizing skin/ stringer debonding failures in reinforced composite panels. Composite Materials: Testing and Design (12th volume), ASTM STP 1274, Deo RB and Saff CR, editors American Society for Testing and Materials; 1996, p. 105–124. [14] O’Brien TK, Chawan AD, Krueger R, Paris IL. Transverse tension fatigue life characterization through flexure testing of composite materials. Int J Fatigue 2002;24:127–45. [15] Lee J, Soutis C. A study on the compressive strength of thick carbon fibre-epoxy laminates. Compos Sci Technol 2007;67(10):2015–26. [16] Lee J, Soutis C. Measuring the notched compressive strength of composite laminates: specimen size effects. Compos Sci Technol 2008;68(12):2359–66. [17] O’Brien TK, Chawan AD, DeMarco K, Paris IL. Influence of specimen preparation and specimen size on composite transverse tensile strength and scatter. NASA report No. ARL-TR-2540; 2001. p. 86. [18] Jiang W-G, Hallett SR, Green BG, Wisnom MR. A concise interface constitutive law for analysis of delamination and splitting in composite materials and its application to scaled notched tensile specimens. Int J Numer Method Eng 2007;69(9):1982–95. [19] Turon A, Dávila CG, Camanho PP, Costa J. An engineering solution for mesh size effects in the simulation of delamination using cohesive zone models. Eng Fract Mech 2007;74(10):1665–82. [20] Hexply 8552 epoxy matrix (180 °C/356 °F curing matrix). Product Data, Hexcel Composites. Publication FTA 072c; October 2008.

F. Hélénon et al. / Composites: Part A 43 (2012) 1017–1027 [21] Camanho PP, Dávila CG, de Moura MF. Numerical simulation of mixed-mode progressive delamination in composite materials. J Compos Mater 2003;37(16):1415–38. [22] Andersons J, König M. Dependence of fracture toughness of composite laminates on interface ply orientations and delamination growth direction. Compos Sci Technol 2004;64(13–14):2139–52. [23] Hélénon F, Wisnom MR, Hallett SR, Allegri G. An approach for dealing with high local stresses in finite element analyses. Composite Part A 2010;41(9):1156–63. [24] Yang Q, Cox B. Cohesive models for damage evolution in laminated composites. Int J Fract 2005;133:107–37.

1027

[25] Harper PW, Hallett SR. Cohesive zone length in numerical simulations of composite delamination. Eng Fract Mech 2008;75(16):4774–92. [26] Hallett SR, Jiang W-G, Wisnom MR. Effect of stacking sequence on open-hole tensile strength of composite laminates. AIAA J 2009;47(7):1692–9. [27] Hallett SR, Green BG, Jiang WG, Wisnom MR. An experimental and numerical investigation into the damage mechanisms in notched composites. Composites Part A 2009;40(5):613–24. [28] Turon A, Camanho PP, Costa J, Renart J. Accurate simulation of delamination growth under mixed-mode loading using cohesive elements: definition of interlaminar strengths and elastic stiffness. Compos Struct 2010;92:1857–64.