Composites Science and Technology 60 (2000) 2849±2856
www.elsevier.com/locate/compscitech
Modelling of splitting and delamination in notched cross-ply laminates Michael R. Wisnom a,*, Fu-Kuo Chang b a University of Bristol, Department of Aerospace Engineering, University Walk, Bristol BS8 1TR, UK Stanford University, Department of Aeronautics and Astronautics, 496 Lomita Mall, Stanford CA94305-4035, USA
b
Received 7 December 1999; received in revised form 12 July 2000; accepted 1 August 2000
Abstract A ®nite-element approach has been developed for modelling the detailed damage development in notched composites. Separate elements are used for each ply, connected together with interface elements to allow delamination between the plies. Interface elements are also used to model splitting at the notch. The approach is applied to a cross-ply laminate with a centre crack loaded in tension, and the results compared with experimental measurements. The model accurately predicts the development of a narrow triangular delamination zone, and the extent of splitting as a function of applied tensile stress. The approach oers scope for improved simulation and understanding of the complex failure processes in notched composites. # 2000 Elsevier Science Ltd. All rights reserved. Keywords: B. Modelling; C. Damage mechanics; C. Delamination; C. Finite-element analysis; C. Notch
1. Introduction Fibre-reinforced composites loaded in tension exhibit notch sensitivity, and this can be an important factor in determining allowable design strains. Considerable research has been undertaken on this subject, and a comprehensive review of the earlier work was presented by Awerbuch and Madhukar [1]. The development of damage is crucial to the behaviour of notched composites, and leads to the well-known hole-size eect, whereby smaller notches have less eect on tensile strength than larger ones. It also gives rise to the phenomenon that residual strength can increase after cyclic loading as a result of the reduction in stress concentration by the fatigue damage. Similar eects are observed whether the notch is a hole or a sharp crack. One approach to this problem is to treat the damage zone as a crack, and apply fracture mechanics. Waddoups, Eisenmann and Kaminski [2] used this method on laminates containing a straight crack by adding an eective ¯aw length to the actual crack length. A good
* Corresponding author.
®t to the experimental data was obtained for tensile strengths of laminates with dierent length cracks. However, the eective ¯aw size is not a material constant, but depends on the layup, as this single parameter cannot represent the dierent forms of damage occurring at the crack tip. Progressive damage modelling allows the nature of the damage to be represented. Chang and coworkers developed an approach including the eects of matrix tensile and compressive failure, ®bre failure, and ®brematrix shearing failure [3,4]. This allowed the extent and type of damage to be modelled as a function of increasing load, and its eect on strength to be predicted. Good results were obtained for the eect of notch size on tensile strength, and the behaviour of dierent layups could be predicted from the same basic material parameters. Delamination was not represented, and in some cases this can have a signi®cant eect. For example in tensile tests on laminates with ®lled holes, certain layups were found to be susceptible to splitting parallel to the edges of the hole, together with delamination [5]. For these layups, bolt clamp-up forces signi®cantly reduced the strength, and this was attributed to suppression of the splitting and delamination that normally alleviate the stress concentration. This phenomenon
0266-3538/00/$ - see front matter # 2000 Elsevier Science Ltd. All rights reserved. PII: S0266-3538(00)00170-6
2850
M.R. Wisnom, F.-K. Chang / Composites Science and Technology 60 (2000) 2849±2856
could not be predicted with progressive damage modelling based only on in-plane failure [6]. The detailed damage arising at notches under tensile loading was studied experimentally by Kortschot and Beaumont [7] and Spearing and Beaumont [8]. For example for a crossply laminate of carbon/epoxy with layup (90/0)s with a centre crack it was shown that splitting occurs in the 0 plies perpendicular to the notch, together with narrow triangular areas of delamination between the 0o and 90o plies, as shown schematically in Fig. 1 [8]. The damage grows with increasing loading until tensile failure occurs. There is also extensive transverse cracking in the 90 plies. A phenomenological approach was developed for modelling the damage. The delamination area was assumed to be triangular with a ®xed angle deduced from the tests, growing in a self similar manner. An energy balance was then used to predict the extent of splitting and delamination under static and fatigue loading, giving very good correlation [9]. An energy approach based on observed damage mechanisms was also used by Johnson and Chang to model delamination from matrix cracks and free edges in angle ply laminates [10]. Analytical models using the critical strain energy release rate criterion were applied to predicting the behaviour of uniformly loaded specimens. Good results were obtained for damage accumulation and strength of angle ply laminates loaded in tension. The eects of ply group thickness, stacking sequence and ply orientation on strength could all be accounted for satisfactorily. The energy balance method can be incorporated into ®nite element analysis for predicting delamination by using the interface modelling approach [11]. Duplicate nodes are placed on either side of any ply interfaces where delamination is expected, and these are connected by elements prescribing the relationship between interfacial stresses and relative displacements between the plies. The approach includes a stress criterion for initiation of damage, and an energy criterion for failure. This means that the analysis determines where failure will occur, whilst allowing propagation to be governed by fracture mechanics. The approach has been used successfully for predicting delamination at free edges
Fig. 1. Schematic illustration of damage in a cross-ply laminate loaded in tension.
[12], cracks [13], discontinuous plies [14] and in fracture toughness tests [15]. In this paper the interface modelling approach is applied to the problem of predicting the detailed damage development in notched composites. Separate planar elements are used for each ply, connected together with interface elements to model delamination between the plies. A similar method is used to model splitting by means of interface elements across the line perpendicular to the notch where the splitting occurs. The approach is applied to a cross-ply laminate with a centre crack loaded in tension, and the results compared with the experimental measurements of Spearing and Beaumont [8]. The development of damage is predicted very well, with good correlation for the extent of splitting and delamination as a function of applied stress. 2. Analysis 2.1. Finite-element model The (90/0)s laminate was modelled with four noded plane stress elements using the ABAQUS ®nite element software [16]. A quarter model was used, with symmetric boundary conditions on both centre planes to represent the eect of the rest of the plate. The width of the model, w/2, was 12 mm, and the length 45 mm. The notch was modelled as a sharp crack with a half width, a, of 4 mm. A relatively ®ne mesh was used adjacent to the notch where splitting and delamination were expected. A ®xed element size with width and length of 0.125 mm and 0.375 mm was used in this area, with larger elements away from the notch. Fig. 2 shows the initial mesh of the upper left quadrant. Since the layup is symmetric, it was only necessary to model half the thickness, i.e. two plies, each taken to be 0.125 mm thick. Separate elements were used to represent the 0 and 90 plies. In the region near the notch with the smallest element size duplicate coincident nodes were used for the two plies to allow relative displacement to occur during delamination. The nodes for the 0 and 90 plies were connected with interface elements, as described in Section 2.2. Away from the notch the same nodes were used for both plies. This approach allows the shear transfer between plies at discontinuities such as ply splits, cracks and free edges to be represented. However, since the elements are plane stress, the eect of the oset between plies in the thickness direction cannot be modelled, and so no mode I loading arises. This is reasonable for the present case, which is dominated by mode II and III shear loading. Duplicate nodes were also used for the 0o plies along the line perpendicular to the edge of the notch where splitting was expected. The 0 elements on either side of the line were attached to dierent nodes, and these were
M.R. Wisnom, F.-K. Chang / Composites Science and Technology 60 (2000) 2849±2856
connected by interface elements to allow the splitting to be modelled. Together with the node for the 90 ply there were therefore three coincident nodes along this line, allowing the splitting and the delamination on either side of the split to be represented. Fig. 3 shows this arrangement schematically. A transverse crack in the 90 ply at the notch can be expected very early due to the high stress concentration and low transverse strength. A crack was therefore assumed to be present right from the start of the analysis. This was done by only applying the symmetric boundary conditions at the plane of the notch to the nodes representing the 0 plies, as shown in Fig. 2. Note that away from the notch where the same nodes were used for both plies, the 90 ply crack was eectively ignored. Orthotropic material properties for T300/914 were de®ned as shown in Table 1 [8]. Initially linear elasticity was assumed, but it became apparent that the non-linear in-plane shear response was very important. This was represented using an equation derived experimentally for a very similar material, XAS/914 [17]:
2851
Table 1 Material properties used in analysis Fibre direction modulus (GPa) Transverse modulus (GPa) Poisson's ratio Shear modulus, linear (GPa) Non-linear shear stress (MPa)ÿ shear strain (%) relation Interface element fracture energy (N/mm) Interface element yield stress (MPa) Interface stiness per unit area (N/mm3)
135.0 9.6 0.31 5.8 =71.41+3.66 ÿ91.52 exp(ÿ ) +20.11 exp(ÿ2 ) 0.4 75 3105
71:41 3:66 ÿ 91:52exp
ÿ 20:11exp
ÿ2
1 Where is the shear stress in MPa, and is the shear strain in %. This is a reasonable assumption since the resin is the same, the volume fraction is 60% in both cases and the materials use similar high strength carbon ®bres. The initial modulus from (1) is 5.5 GPa, close to the assumed linear value of 5.8 GPa for T300/914 [8]. The non-linear response was implemented in ABAQUS
Fig. 2. Finite element mesh and boundary conditions.
Fig. 3. Use of duplicated nodes and interface elements.
2852
M.R. Wisnom, F.-K. Chang / Composites Science and Technology 60 (2000) 2849±2856
by means of a UMAT user subroutine based on elasticity, which is satisfactory for monotonic loading. The shear response was assumed to be independent of the ®bre direction and transverse stresses. The model was suciently long to reach a uniform tensile stress away from the notch. All nodes at the end of the specimen were therefore constrained to have the same displacement in the loading direction, and the load was applied by means of a single force at one node. 2.2. Interface modelling The interface was modelled using non-linear springs between the coincident nodes representing the 0 and 90 plies, with the characteristics shown in Fig. 4. The springs are initially assumed to be elastic, with a stiness based on the characteristics of the resin rich layer between the plies. A spring stiness per unit area, k, of 3105 N/mm3 was used, derived from typical values of shear modulus for epoxy of 1.5 GPa and resin layer thickness of 0.005 mm. When the stress reaches a critical value, it is assumed that relative deformation can take place between the plies at constant stress, representing an elastic-perfectly plastic interface response. A critical stress of 75 MPa was used, based on a typical epoxy yield stress. The sensitivity of the results to these parameters is discussed in Section 3.4. When a critical relative displacement is reached, the interface is assumed to fail. The area under the shear stress- relative displacement graph corresponds to the delamination fracture energy. A value of 0.4 N/mm was used based on previous work [9]. If relative motion could only take place in one direction, a single degree of freedom non-linear spring element could have been used to represent this behaviour. However, in the present case relative deformation can occur in both the x and y directions, and these need to be coupled together. This was achieved by means of a two noded user de®ned element in ABAQUS with both x and y degrees of freedom at each node. The two interlaminar shear stresses xz and yz are calculated from the relative displacements in the x and y directions multiplied by the spring stiness per unit area, and hence the resultant shear stress is determined. If this is
below the yield stress the interface is assumed to be elastic, with stinesses in both directions equal to the interface stiness per unit area multiplied by the eective area at that node. In this study all the elements were rectangular with four nodes. The interface area for each node was therefore taken as the sum of one quarter of the area of every element attached to that node. After the yield stress has been reached, it is assumed that relative displacement can occur with no change in the resultant shear force. The stinesses of the springs in both the x and y directions are set to zero. The direction of relative displacement is allowed to vary, with the assumption that in each increment the forces in the x and y directions are proportional to the x and y displacement increments. The magnitudes of the resultant displacements in each increment are summed and when the critical value is exceeded, the interface is assumed to fail, and the forces are set to zero. The critical displacement is taken as the value of the interfacial fracture energy divided by the yield stress. When delamination occurs there is also energy released due to relaxation of thermal residual stresses. However, calculations showed this to be small in comparison with the fracture energy, and so residual stresses were neglected in the analysis. For convenience the same interface element was also used to model splitting in the 0 ply. In this case there is really only one degree of freedom at each node: the shear displacement tangential to the split. However, the displacements perpendicular to the split are constrained by the sti ®bres in the 90 ply, and so no signi®cant relative motion can occur in this direction prior to delamination. The interface elements with degrees of freedom in both directions therefore behave in a very similar way to single degree of freedom springs, since splitting occurs before delamination. The area associated with each node is now the 0 ply thickness multiplied by the eective length. This was taken as the sum of half the lengths along the split of the two elements attached to that node. The last node at the notch only attaches to one element, and so the area for this interface element was halved. Since the splitting is dominated by mode II shear deformation, the same yield stress and fracture energy parameters were used as for the interface elements modelling the delamination, which is also shear dominated. 2.3. Solution
Fig. 4. Assumed interface characteristics.
A load equivalent to a stress of 500 MPa averaged over the complete width and thickness was applied in increments. An iterative solution was required at each increment due to the non-linear shear material response and the non-linearity introduced by the interface elements. In most cases the solution converged within six iterations. However, at certain points convergence
M.R. Wisnom, F.-K. Chang / Composites Science and Technology 60 (2000) 2849±2856
became more dicult, requiring more iterations or a smaller increment size. This was aggravated by the discontinuity arising when an interface element failed, at which point the corresponding loads were suddenly removed. The ABAQUS option for automatic solution control allowed the increment size to be reduced where necessary, giving convergence in most cases. Occasional diculties were still experienced at the point where an element was about to fail, with the solution switching between the element in the intact and failed states on successive iterations. In such cases the oending element was removed to enable the solution to be taken further. The length of the split and the extent of the delamination were determined from the number of interface elements that had already completely failed or were just about to fail. Since failure occurs in discrete steps of one element at a time, the accuracy of this procedure is limited by the element size. 3. Results and discussion 3.1. Initial results Splitting initiated from the end of the notch at a relatively low load, and then progressively increased in length. This was accompanied by the formation of a narrow zone of delamination approximately triangular in shape, as observed experimentally [8]. Fig. 5 shows the split length as a function of the average stress applied to the specimen compared with the experimental results, and the correlation is good. Fig. 6 shows the exaggerated deformed mesh at an applied load of 367 MPa, and indicates the pattern of delamination. Away from the notch the 90 and 0 ply elements are coincident, but where there is damage the dierent elements can move relative to each other, and the mismatch of element boundaries shows the extent of the delamination. Similarly the splitting can be seen from the mismatch in elements along the line perpendi-
Fig. 5. Correlation of predicted and measured splitting.
2853
Fig. 6. Exaggerated distorted geometry at 367 MPa applied stress.
cular to the notch. The width of the delaminated region at the notch is 0.625 mm, and the length is 9.56 mm, giving an average angle of 3.7 . This compares closely with the angle of about 4 for the triangular delamination measured experimentally [8]. The extent of delamination is slightly less than the split length, by one node inside the notch (0.375 mm), and by three nodes outside (1.125 mm), showing that it is the splitting that drives the delamination. The shape of the delaminated region is similar at lower loads, indicating self similar damage growth, as seen in the tests [8]. The yield zone extends two nodes beyond the tip of the split, i.e. 0.75 mm. At the notch the yield zone extends across the width beyond the delaminated area by one node inside the notch and by three nodes outside, i.e. a total additional width of 0.5 mm. The yield zone has just reached the edge of the region with duplicated nodes, and so the re®ned mesh would need to be extended in order to take the analysis further. 3.2. Eect of mesh size The analysis was repeated with the element size doubled in both directions in order to investigate the sensitivity of the results to mesh re®nement. The element size near the notch was now 0.75 mm long and 0.25 mm wide, compared with 0.3750.125 mm in the original model. The interface element areas were adjusted accordingly. All other aspects of the analysis were kept the same. Fig. 7 shows the results for the predicted split length for both meshes, and they are similar. The coarser mesh gives a slightly less smooth response, since there are fewer nodes to represent the growing damage. The biggest discrepancy occurs at the beginning, where the coarse mesh does not have sucient nodes to model the sharp angle of the delamination pattern. At higher loads the triangular delamination pattern is very similar to that obtained before. The predicted angle is 3.3 at
2854
M.R. Wisnom, F.-K. Chang / Composites Science and Technology 60 (2000) 2849±2856
327 MPa applied stress, close to the value of 3.7 with the ®ner mesh. These results show that the analysis is not sensitive to the mesh size provided there are sucient elements to represent the damage. This is in contrast to models using stress based criteria for damage, which can show signi®cant dierences with dierent element sizes. 3.3. Eect of material modelling The analysis was repeated using a constant shear modulus of 5.8 GPa in order to investigate the eect of shear non-linearity. The delamination pattern was much broader, and so it was necessary to use a dierent mesh. The revised mesh had the same number of elements across the width, but with an element width of 0.25 mm rather than 0.125 mm near the notch. The mesh size along the length was 0.75 mm, giving the same aspect ratio as before. The dierent mesh used for this case is not expected to in¯uence the results greatly given the insensitivity to mesh size shown in Fig. 7. The predicted split length is compared with the original baseline analysis in Fig. 8. The results are quite dierent,
Fig. 7. Eect of mesh re®nement on predicted splitting.
Fig. 8. Eect of material properties on predicted splitting.
with much less splitting for the case with linear shear response. The delaminated area has a similar triangular shape to before, but the angle is now about 15 , much larger than the 3.7 obtained previously. These results do not correlate with the experimental measurements, and show the importance of taking account of the shear non-linearity in the analysis. This is not surprising, as shear stresses in excess of 250 MPa were obtained in the 90 ply with the linear analysis, well above what the material is capable of carrying. Transverse cracking also takes place in the 90 ply [8], reducing the eective modulus. To investigate the sensitivity of the results to this eect the analysis was repeated with the transverse modulus halved in both the 90 and 0 plies. Shear non-linearity was included, and all other aspects were the same as for the baseline case. Fig. 8 shows the results, and it can be seen that they are very similar to the original ones. The splitting and delamination behaviour should therefore not be sensitive to the eect of transverse ply cracking, and it is reasonable to neglect it. 3.4. Eect of interface element parameters The initial elastic stiness of the interface elements was based on typical characteristics of the resin layer between plies. To investigate the sensitivity of the results to the value used, the analysis was repeated with the spring stiness halved. Everything else was kept the same. Fig. 9 shows that the results overlie the original ones, demonstrating that the precise value chosen for this parameter is not critical. The analysis was also repeated with the value of the yield stress for the interface reduced from 75 MPa to 37.5 MPa. These results are also shown in Fig. 9, and again, the progression of damage is similar to the original case. The size of the yield zone has greatly increased, extending 4.5 mm beyond the end of the split at 300 MPa applied load. There is also an increase in the width
Fig. 9. Eect of interface element parameters on predicted splitting.
M.R. Wisnom, F.-K. Chang / Composites Science and Technology 60 (2000) 2849±2856
Fig. 10. Eect of splitting fracture energy.
of the yield zone, which at this load extends beyond the delaminated area by a total of 0.75 mm at the notch. Changing the fracture energy would greatly aect the results, as this is the key parameter which controls the damage. However, it can be seen that the results are not very sensitive to the values of the other two interface element parameters. 3.5. Eect of splitting fracture energy In previous analysis of these notched tensile tests a value of 0.158 N/mm was used for the splitting energy, corresponding to mode I rather than mode II fracture [9]. In this study the same value of 0.4 N/mm was used for the splitting and delamination energy since it was considered that both damage mechanisms should be controlled by shear fracture. To investigate the eect of this assumption, the analysis was repeated with the splitting fracture energy reduced to 0.158 N/mm. Fig. 10 shows the results. In the early stages of damage progression the reduced splitting energy gives an improvement, matching the experimental results very closely. For longer split lengths, the original analysis gives better correlation. These results suggest that initially mode I may make a greater contribution to the splitting. However, the ®nite element analysis indicates that the transverse stresses in the 0 ply are much lower than the shear stresses, even in the early stages of damage development, supporting the assumption of mode II dominated fracture. Residual thermal stresses may have some in¯uence on this behaviour, and further investigation is required. 4. Conclusions An approach has been developed for predicting the damage development in notched composites. The model
2855
represents the physical processes of splitting and delamination, and is therefore able to simulate closely the behaviour observed experimentally. Results for the split length as a function of applied tensile stress correlate very well with experimental measurements on notched crossply laminates. The self similar damage growth and triangular delamination zone with an angle of about 4 also match the experimental observations. The predicted damage is largely controlled by the mode II fracture energy, and is not sensitive to the mesh size, the assumed interface stiness and yield stress or the eect of transverse ply cracking. Nonlinear shear behaviour is important, and must be included in the analysis. The approach allows the eect of damage in redistributing the loads around stress concentrations to be modelled accurately. It can be applied to other geometries and layups, and oers scope for improved simulation and understanding of the complex failure processes in notched composites. Acknowledgements This work was supported by the UK Engineering and Physical Sciences Research Council under grant GR/ M90214.
References [1] Awerbuch J, Madhukar MS. Notched strength of composite laminates: predictions and experiments Ð a review. Journal of Reinforced Plastics and Composites 1985;4:3±159. [2] Waddoups ME, Eisenmann JR, Kaminski BE. Macroscopic fracture mechanics of advanced composite materials. Journal of Composite Materials 1971;5:446±54. [3] Chang FK, Chang KY. A progressive damage model for laminated composites containing stress concentrations. Journal of Composite Materials 1987;21:834±55. [4] Chang KY, Liu S, Chang FK. Damage tolerance of laminated composites containing an open hole and subjected to tensile loadings. Journal of Composite Materials 1991;25:274±301. [5] Yan Y, Wen WD, Chang FK, Shyprykevich P. Experimental study on clamping eects on the tensile strength of composite plates with a bolt-®lled hole. Composites 1999;30A:1215±29. [6] Chang FK, Qing X, Sun H-T, Yan YM. Damage tolerance-based design of bolted composite joints, Final report to Boeing Commercial Airplane Group under contract No. TZ-370444-07LHN and to FAA Technical Center under contract No. DEPT/95-G012-0001, Stanford University, Department of Aeronautics and Astronautics, 1999. [7] Kortschot MT, Beaumont PWR. Damage mechanics of composite materials. I: Measurements of damage and strength. Composites Science and Technology 1990;39:289±301. [8] Spearing SM, Beaumont PWR. Fatigue damage mechanics of composite materials. I: Experimental measurement of damage and post-fatigue properties. Composites Science and Technology 1992;44:159±68. [9] Spearing SM, Beaumont PWR, Ashby MF. Fatigue damage mechanics of composite materials. II: A damage growth model. Composites Science and Technology 1992;44:169±77.
2856
M.R. Wisnom, F.-K. Chang / Composites Science and Technology 60 (2000) 2849±2856
[10] Johnson P, Chang FK. Characterization of matrix crack-induced laminate failure, part II: analysis and veri®cations, Journal of Composite Materials, in press, 2000. [11] Cui W, Wisnom MR. A combined stress based and fracture mechanics based model for predicting delamination in composites. Composites 1993;24:467±74. [12] Schellekens JCJ, de Borst R. A nonlinear ®nite-element approach for the analysis of mode-I free edge delamination in composites. Int J Solids and Structures 1993;30:1239±53. [13] Wisnom MR. Modelling the eect of cracks on interlaminar shear strength. Composites 1996;27A:17±24.
[14] Petrossian Z, Wisnom MR. Prediction of delamination initiation and growth from discontinuous plies using interface elements. Composites 1998;29A:503±15. [15] Mi Y, Cris®eld MA, Davies GAO, Hellweg HB. Progressive delamination using interface elements. Journal of Composite Materials 1998;32:1246±72. [16] ABAQUS version 5.7-3, Hibbitt, Karlsson and Sorensen Inc., 1080 Main Street, Pawtucket, RI 02860, USA. [17] Wisnom MR, Haeberle JG. Prediction of buckling and failure of unidirectional carbon ®bre-epoxy struts. Composite Structures 1994;28:229±39.