- Email: [email protected]
Shear driven delamination propagation in two dimensions
PII: S1359-835X(97)00015-8
Composites Part A 28A (1997) 757–765 q 1997 Published by Elsevier Science Limited Printed in Great Britain. All rights reserved 1359-835X/97/$17.00
Shear driven delamination propagation in two dimensions
G. A. O. Davies*, P. Robinson, J. Robson and D. Eady Department of Aeronautics, Imperial College of Science, Technology and Medicine, Prince Consort Road, London SW7 2BY, UK (Received 13 January 1997; revised 20 January 1997)
This paper investigates, theoretically and experimentally, the propagation of delamination in shear mode II between internal ply surfaces in a laminated composite plate. By considering two special cases the problem is approximated as an axisymmetric model amenable to analytical solutions. In one case the circular crack front expands and in the other it contracts. It is shown clearly that only a (linear) fracture mechanics approach will explain the mode and threshold of failure, but that a strength-based criterion can be used to initiate propagation without the need for an initial flaw. Apart from explaining some curious phenomena, the simple axisymmetric analytical solutions may be used as benchmarks for testing the ability of finite element models to predict plate delamination when the crack front is curved (at present only simple beam models are available for the standard DCB and ENF tests, but even these tests do not really develop straight crack fronts). q 1997 Published by Elsevier Science Limited. (Keywords: delamination; mode II; finite element simulation)
INTRODUCTION Laminated carbon-fibre epoxy composites have been a success story whenever high specific strength and stiffness is needed. Military and civil aircraft, and racing cars are clear examples, although the total yearly production of carbon composite structures owes far more to sporting goods than transport vehicles. Nevertheless, the use of composites is still on the increase in aerospace for example, where the manufacturing costs are coming down by a combination of automated production methods and far fewer component parts compared with metal assemblies. It has been estimated that a new generation civil aircraft wing of carbon composites can now be made which is 20% cheaper, 20% lighter, and have an infinite fatigue life compared with the metal equivalent. The main weakness of laminated composites is still the very low through-thickness strength which can be less than one fortieth of the in-plane values. Furthermore both the fibres and resin in all composite systems are brittle when compared to ductile metals, so even small stress concentrations can initiate a very local failure which may then propagate in an unstable fashion. Obvious weaknesses and potential failure sources will be between bonded components or at any form of discontinuity, in load or geometry, but another hidden menace is the interior delamination due to shear stresses which will likely be a maximum in the mid-thickness region. Such * To whom correspondence should be addressed
internal delamination may not be serious for laminates loaded by in-plane tension but they are potentially disastrous for plates and shells designed to resist compression. It is for such structures, where the specific stiffness is the design criterion, that carbon composites come into their own. The history of full-scale composite structure testing has not been altogether a happy one, and a significant number of aerospace structures have failed to meet their design load due to premature through-thickness induced failure. Such weaknesses have to be overcome by either improved resin properties, through-thickness reinforcement, or more importantly by better detailed design. It is therefore vital that designers have the analytical tools available for predicting the debilitating effect of through-thickness stresses and for evaluating them in the first place. Commercial finite element codes have become very user-friendly, are equipped with better composite models, and the modern workstation is giving practitioners the capability of modelling highly complex stress fields is in 3-D regions. It is important therefore to know what to do with the answers. This paper addresses just one source of failure driven predominantly by interlaminar shear, and the commonest situation in plates where the delamination front is curved to enclose a 2-D delaminated area.
DELAMINATION FAILURE CRITERIA It is well known that using simple strength criteria, possibly couched in terms of strains, can be misleading for brittle
757
Shear driven delamination propagation in two dimensions: G. A. O. Davies et al. materials. This is particularly the case in the presence of rapidly varying stress fields. A fracture mechanics or a damage mechanics approach may therefore be a better model. However, strength criteria have traditionally been the most widely used, and the most common test for interlaminar strength is still the 3-point bending test using a short beam1 which does give reasonably consistent values for the interlaminar shear strength. One reason for its success is the fact that the rig is usually soft loading and the shear force is constant along the beam; so that once a failure initiates, the delamination is virtually instant. Even so, the initiation of failure in the complex stress field, local to the load and supports, is a challenge to model2. Any failure in a brittle material has to capture the physics properly. Stress concentrations will not be relieved by plastic flow and yet stress concentrations can be exceptionally high in laminated composites, even in the absence of clear stress raisers—for example edge effects. Moreover, strength-based criteria alone will not scale; geometrically similar structures of different size and with identical stress fields will fail at different stress thresholds. Fracture mechanics will often resolve this dilemma since it does directly address the physics of crack propagation and the balance between the structure’s strain energy release and the energy absorbed by the movement of the advancing crack front. The nature of the stress field into which that crack propagates will also need to be correctly modelled. However, fracture mechanics needs the insertion of a flaw before the prediction of crack propagation can take place and the energy release rate (or stress intensity factor) can be quite sensitive to the size of the flaw length l. Classical linear elastic fracture mechanics in isotropic materials shows strength to be proportional to l ¹1/2. If we assume that all composite structures contain natural flaws, and then use a fracture criteria to explain experimentally found strengths, we find the initial flaw to be of the order of millimetres,3 whereas microscopic examination fails to reveal flaws of such sizes4. Nevertheless, we need separately to emulate initiation and propagation (stable or unstable) consequently and the strength criteria would seem logical before the crack appears, and a fracture criterion thereafter. This strategy has been suggested by many workers5–8 and seems to work reasonably well, the main problem being the expense and the monitoring of the crack front. By far the most common approach to compare numerical models and experimental results has been to choose 1-D tests such as the double cantilever beam (mode I), the end notch flexure test (mode II) and the short beam under 3 or 4-point bending. Even these supposedly simple cases are not strictly 1-D and the crack front can develop a curved profile at the plate free edges9,10. The aim of the fracture approach is to evaluate numerically or analytically the energy release rate G when the crack front is allowed to move a small amount, and then compare this with the known critical value. The critical value may be for a pure opening mode I (G IC), for a pure shear mode II (G IIC) or for a mixed mode11. The usual way of finding G is by virtual crack closure12 or virtual crack extension13 and hope that small crack perturbations are self similar. It appears that quite small front perturbations, using
758
ultra-fine or adapting finite element meshes, can be needed, although several workers10 have successfully propagated curved and developing crack fronts along which the perturbations are chosen to ensure propagation only when G rises to its critical value. A completely different approach, which avoids tracking a crack front completely, is the use of interface elements14 in which the interface stiffness has a loading portion which rises to a ‘strength’ and decays to zero so that the area under the load–displacement curve is equal to G C. This strategy therefore models initiation and propagation and does not have a true discrete front advancing through the mesh. Unfortunately, iterative solutions of structures with softening can lead to numerical problems which need a fair amount of computing effort to overcome11. However it is not the purpose here to debate these strategies but to validate the initiation-propagation strategy for 2-D problems, where both the energy release rate and the crack front curvature are changing simultaneously.
1-D DELAMINATION IN BEAMS The most tested benchmark cases are the standard tests for finding mode I and II fracture toughness. We concentrate on the shear mode II, or Edge Notch Flexure (ENF) and End Load Split (ELS) tests15 (see later in Figure 10). In the latter case if the two arms of the split cantilever are identically balanced and the equivalent flexural modulus of the chosen laminate is E, then the applied tip load P is shared equally and produces a pure shear field at the end of the crack length a. Simple beam theory gives the tip deflection in terms of a 3 and the energy release rate, after allowing a perturbation da, as: G¼
9P2 a2 4B2 Eh2
(1)
where B is the beam width and h the depth of each split half. This energy is solely due to bending stresses (although it can be extended to include finite rotation16) and the bending stresses increase linearly away from the applied load so that the crack is advancing into an increasing stress field. The energy release therefore increases with a in equation (1) and the delamination is therefore unstable under the constant dead load P. This test is an accepted standard but the theoretical model is 1-D, and the crack front does not change in length. We now turn to two examples where this is not so.
2-D DELAMINATION IN THIN PLATES To keep the analysis simple, and free from further complications in numerical strategies, we pose two axisymmetric problems with very different internal shear fields. There is of course no such structure as a laminated composite plate of constant thickness and having axisymmetrical properties, so this is our first approximation. However, the tests described will all deploy quasi-isotropic lay-ups with a large number of plies up to 32 for thicker plates so the structural response will be close to axisymmetrical. The first case is
Shear driven delamination propagation in two dimensions: G. A. O. Davies et al. shear strength (t u) then suppose the laminate has a large number of plies in the quasi-isotropic lay-up, none of them blocked together. The maximum shear stress will then be approximately 1.5 3 the mean, and we can predict failure when: 1:5
P ¼ tu 2prt
and the circular delamination area will then be given by A ¼ pr2 ¼
Figure 1 Variations in shear and peeling stresses away from a point load
that of a single concentrated force where the far-field through-thickness stress decays like 1/r and the second case a uniform pressure, reacted by a far-field boundary, where the shear stress increases like r. Concentrated force It will be helpful to consider the nature of the mean internal stresses j zz and t rz as shown in Figure 1. The farfield shear stress t rz is simply given by P/2prt until we approach the axis of symmetry under the load where, by symmetry, it must be zero. Clearly the stress and strain field immediately under the contact force is highly complex and nonlinear, but the stresses will approach the far field values in a transition region of order t as indicated in Figure 1. From the equilibrium equations it is found that the throughthickness peel stress j zz is proportional to the derivative of the shear and is sketched in the figure with a compressive zone under the applied force and a tensile zone further away, in which numerical models have shown the stress to be small in magnitude. The zone vulnerable to delamination is therefore away from the load when the compressive stresses change to tensile and the shear is still high. Highly detailed finite element models confirm this ‘average’ description and indeed tests show sometimes a delaminated annulus surrounding an intact centre17. If we wish to use a strength criteria, based on interlaminar
2
9 P 2 t 16pt u
(2)
Because the shear stresses decrease to zero in the vicinity of the origin, we do not expect equation (2) to be valid for a distance of, say, r ¼ t. A somewhat arbitrary bound on this model is therefore noted as A ¼ pt 2. As it happens a very large number of low velocity impact tests, on an equally large variety of plate sizes and thicknesses, were available for testing equation (2). The material was HTA/6376 for all plates. These tests have been summarised in Davies and Zhang17 but for completion Figure 2 has been reproduced here to see if equation (2) is a valid model. In low velocity impact, the force generated is many times greater than the distributed plate inertia forces, so the above static model is fair. However, we admit that delamination is far from a single circular patch. The damage is not confirmed to a single plane—in fact the middle plane of a symmetrical stacking sequence is not even a true interface. Instead the delamination consist of elongated oblongs in the direction of the ply fibres at the delamination lower surface, and these delamination may occur on four or more interfaces in the central region of the plate thickness. However the envelop, as revealed in ultrasonic scans, is nearly a perfect circle17. It appears that the envelop of these oblong fronts does advance, as a circle for quasi-isotropic laminates, with increasing load. With these reservations Figure 2 shows the application of equation (2) to a series of impacted plates of many sizes and having 8, 16, and 32 plies. An interlaminar shear strength of 50 MPa for HTA/6376 composite is a lower bound to take into account the fact that delamination occurs between different ply orientations. We see that the three strength curves, even with the lower bound imposed, do not reflect the behaviour of the damage. For the thinnest plates of 1 mm the map is reasonably up to a force of 1000 N, but for the thickest plate equation (2) completely misses the very rapid development of damage which occurs at an impact force of approximately 5500 N. In fact it was not possible in these tests to control the damage near this load, the development of delamination was truly unstable, and suggested therefore a fracture mechanism. If we therefore turn to Figure 3, taken from Davies et al.18, we can evaluate energy release rate, around the boundary of a single circular delamination, in the same way as equation (1) using isotropic plate theory for the deformation. If the delamination did not exist in the plate zone bounded by A and A9B9 we can easily find the deflection of the load P due solely to the strains in the region ABB9A9 and taking the normals like AA9 and BB9 as fixed in space. Denote this
759
Shear driven delamination propagation in two dimensions: G. A. O. Davies et al.
Figure 2
Strength-based predictions for delamination area. (X denotes lower bound)
therefore the true extra deflection of P due to the insertion of the delamination, and the net work available to drive the crack is given by: 1 9Pa2 (1 ¹ v2 ) U ¼ P(du ¹ dd ) ¼ 2 64pEh3 Notice that this energy is a function of a 2, as for the 1-D beam equation (1), but the crack front length increases with a, so that when we equate: G·2pa da ¼
]U da ]a
we find G¼
Figure 3 Simplified model of internal delamination for a quasi-isotropic
If there is a fraction-free delamination of radius a then the two separate halves move together and share the load P equally, whence the ‘damaged’ deflection is this time given by: P dd ¼ 3 a2 (1 ¹ v2 )4ph3 2 The normals AA9 do of course rotate, but the rotation is due to the strains in the undamaged region outside ABB9A9, and is the same for both cases. The difference d d ¹ d u is
760
(3)
Equation (3) is remarkable in that it is independent of the size of the delamination, a. (In fact, in view of the difficulties in verifying crack front positions in tests, it might be a candidate for replacing the ELS test!). It does successfully explain the curious behaviour of the damage maps in Figure 2, if we equate equation (3) to the critical value of G IIC and invert. The critical force level
‘undamaged’ deflection as: du ¼ 3Pa2 (1 ¹ v2 )=4p(2h)3
9P2 (1 ¹ v2 ) 64p2 Eh3
P2c ¼
8p2 E(2h)3 GIIC 9(1 ¹ v2 )
(4)
and is also independent of crack size. Thus there is a genuine threshold force, below which no damage will propagate, but when the threshold is reached the amount of damage is indeterminate if the threshold force is kept constant. There is nothing to be gained by trying to use fracture analysis to evaluate the extent of damage. It also explains why there is so much scatter in damage maps. It is not due to variable non-reproducible composite fibre and resin properties, as is often suggested.
Shear driven delamination propagation in two dimensions: G. A. O. Davies et al.
Figure 4
Fracture-based predictions for force threshold
In Davies and Zhang17, equation (4) was used with an average value of G IIC of 0.8 N/mm to account for the differing interfaces at which delamination occurs. Figure 4 shows that the thresholds are predicted reasonably well. For the set of thickest two plates the strength-based initiation is not activated at all, and the lower bound from Fig. 3 is believable. The unstable fracture threshold is however predicted accurately. For the thinnest set of plates, the lower bound does not now inhibit the strength-based curve as much and the tests do show the parabolic rising curve until an unstable growth at about 1000 N. The fracture-predicted threshold value on this occasion is a little on the low side.
The point to note is that, although the shear stress decreases away from the outer boundary, so does the perimeter of an advancing crack front. It is shown19 that as the crack advances from the outer edge of the circular plate of radius R, by a radial distance l, the mode II energy release rate is GII ¼
9 p2 (1 ¹ v2 ) [(2R ¹ l)l]2 32 Et3
Inverting this equation there should be a critical pressure, given in terms of G IIC, as
4 2GIIC Et3 1 Pc ¼ 2 3 (1 ¹ v ) (2R ¹ l)l
Uniform pressure The other case of a plate loaded by uniform pressure, and needing therefore an outer edge support, has completely different stress field, again amenable to simple theory. If the applied normal pressure is p, the average shear stress at any radius distance r from the centre, will be given by: ppr2 =2prt and allowing for the maximum stress t to be 1.5 3 the mean, we have: t¼
3 pr 4 t
(5)
Equation (5) can be used to estimate the radius of the inner boundary of delamination by putting t equal to the interlaminar shear strength. As the shear increases linearly the delamination will extend from the outer support to the front but will apparently be a stable delamination since the stress decreases from the delamination front into the undamaged central region. The true test of stability needs an energy release approach. In this case the work is done by a distributed pressure, but it is amenable to an energy release analysis similar to the previous example19.
(6)
(7)
Unlike the previous example this pressure loading does depend on the delamination length l and for small l p R, it varies inversely like 1/l. As l increases the critical pressure in equation (7) decreases, indicating that the delamination is unstable, and not as the strength criteria would have us believe. This is not too dissimilar from in-plane classical fracture phenomena where the applied critical tension varies like (crack length) ¹1/2, so we have returned to the dilemma of selection an initial flaw. We therefore have conducted a series of experimental tests to investigate the true nature of the failure mode, initiation, and development of the crack fronts.
EXPERIMENTAL DETAILS To pick up possible scale effects, two sizes of circular plates were chosen, 130 m diameter and 270 mm diameter. The plates, this time made of T800/924 carbon epoxy, were held with the edges clamped and an artificial delamination using 0.0125 mm thick Teflon insert was placed in the mid-plane during fabrication. In a symmetrical lay-up the mid-plane is not a true interface, so to represent a natural delamination
761
Shear driven delamination propagation in two dimensions: G. A. O. Davies et al.
Figure 5
Figure 6
Test rig for applying uniform pressure
History of displacement and pressure in a typical test
source the adjacent plies were chosen to be 6 45 by using the following 32-ply (nominal 4 mm) lay-up [A= ¹ A=A= ¹ A] where A ¼ 45/90/ ¹ 45/0/0/ ¹ 45/90/45.
762
This lay-up is balanced with no twist–bend coupling and was perfectly flat when removed from the autoclave. Delamination widths of 5, 10, 15, 20, and 25 mm were chosen. The plates were clamped by bolts and gasket as shown in Figure 5, which sealed the container of water to be pressurised using a 200 bar hand-pump. The water pressure and the plate’s centre deflection were monitored continuously during the manual loading. The loading was increased at as near constant rate as possible until the delamination was clearly heard to propagate. ‘Failure’ was instant and unambiguous, and the pressure dropped immediately with very little change in plate deflection since water is virtually incompressible. Figure 6 shows a typical history. After failure, three post-mortem diagnoses were used. Each plate was ultrasonically C-scanned to confirm the unstable fracture, and Figure 7 shows a typical output. The original is more clearly colour coded, to reflect the time-offlight to the delamination depths, and shows that the artificial (5 mm) delamination has propagated into the plate centre over the regions subtending roughly 10 to 2 o’clock
Shear driven delamination propagation in two dimensions: G. A. O. Davies et al.
Figure 7 Final delaminated zones starting from 5 mm around the edge. The grey scale is the mid-plane
and 3 to 9 o’clock. Away from the axes of symmetry the 6 45 interface becomes 0/90 to an advancing circular front, and here the delamination switched to the weaker 90/45 interface. Nevertheless, the unstable front is almost axisymmetrical and probably would have advanced to the place centre had the pressure loading between sustainable. Fractographic examination of the failure surfaces showed characteristic cusps indicating a dominant mode II fracture mostly in the radial direction. A very small region of mode I peel was visible at the edge of the Teflon insert. Finally it was predicted that bending stresses around the clamped edges would exceed 2000 MPa for the large plates. On taking sections the polished micrographs confirmed extensive compressive microbuckling near the top surface and cracking near the lower (Figure 8). The consequent degradation of the plate’s bending stiffness would modify the clamped boundary condition, but fortunately this does not affect the simple theoretical predictions of equations (5) and (7). Turning to these equations we need to select an interlaminar shear strength, a fracture toughness, and a flexural modulus. Unfortunately, all three vary with orientation around the plate. There are also different values of fracture toughness for initiation and propagation. We therefore selected upper and lower bounds of: Interlaminar Shear Strength: 120 and 80 MPa G IIC: 1.1 and 0.55 N/mm E: 60 and 44 GPa These bounds are rather widely spaced, although as indicated in equation (7) the critical pressure depends on the square root of the product of toughness and modulus. Nevertheless the results plotted in Figure 9 do appear conclusive for both plate sizes. The shear-strength based criteria does predict the initiation of failure for the undamaged plates, but becomes progressively inaccurate as the initial delamination size increases. The fracture-based
Figure 8 Micrograph near top and bottom surfaces under high bending strains
criterion is quite respectable in spite of the simplifying nature of our theory. Thus the experiments do confirm the nature, mode, and critical pressure of the theoretical failure model. The test data clearly switches from strength to fracture even though the transition is not quite sharp as the two separate theories would indicate. In this case, rather like the 1-D short beam test, the shear strain does not vary drastically for a small distance from the plate edge, so there is no need to explain why an undamaged plate needs an initial flaw for the fracture criteria to be triggered. It is permissible to use the strength predictor until the delamination length causes the fracture criteria to override it. CONCLUSIONS • Without introducing the complexities of a numerical analysis it has been possible to select two (quasi) axisymmetrical problems of plate bending where the failure mode is delamination, driven by mode II shear, and the crack front envelope is nearly circular. • In both cases a strength criterion can be used to initiate a delamination, but the energy release model is needed to
763
Shear driven delamination propagation in two dimensions: G. A. O. Davies et al.
Figure 9
(a) Small plate failure. (b) Large plate failure
predict the subsequent growth and the stability of that growth. • A numerical simulation using interface elements should be successful if the load/deflection law increases to a maximum strength but thereafter decreases so that there is total energy release equal to G IIC. • A fracture based virtual crack closure (or extension) strategy should simulate the advancing crack front, but it is suggested that the initial node releasing may be based on interlaminar strength until (and if) the fracture analysis indicates growth at a lower load than the strength based criterion. • Finally we note that there is much current research aimed at using Finite Element models to initiate and propagate delamination. In reality most crack fronts are curved,
764
even for the simple beam tests illustrated in Figure 10, and the F.E. community has to validate the numerical models with experimental observation involving all the usual uncertainties of strength, modules, toughness, and the position of the crack front. We therefore offer these two axisymmetrical solutions as a benchmark for curved crack front propagation.
ACKNOWLEDGEMENTS The authors would like to thank the DRA (Structural Materials Centre, Farnborough) for supporting projects on low velocity and ballistic impact which revealed the need for these basic failure studies.
Shear driven delamination propagation in two dimensions: G. A. O. Davies et al. 6. 7. 8. 9. 10.
11. 12.
Figure 10 Standard ENF and ELS tests for mode II fracture toughness
13. 14.
REFERENCES 1.
2. 3. 4. 5.
Curtis, P.T., ed., Crag test methods for the measurement of the engineering properties of fibre-reinforced plastics. Royal Aerospace Establishment, Tech. Rep. 88012, Farnborough, UK, February 1988. Wisnom, M.R., Cui, W.C. and Jone, M., Failure mechanisms in 3 and 4-point short beam bending tests of u-d glass/epoxy. Journal of Strain Analysis, 1992, 27(4), 235–243. Davies, G.A.O. and Cui de, Yu., Three dimensional edge effects in composite plates. Computers and Structures, 1980, 31(1), 11–16. Hull, D., An Introduction to Composite Materials. Cambridge University Press, 1981. Bostaph, G.M. and Elber, W., A fracture mechanics analysis for delamination growth during impact on composite plates. ASME Winter Conference, 1982.
15. 16. 17. 18. 19.
Wisnom, M.R., Interface element approach to modelling mode I and II fracture. In ICCM, 10, Whistler, Canada, 1995, pp. I:77–83. Minguet, P.J. and O’Brien T.K., Analysis of composite skin/stringer bond failure using a strain energy release rate approach. In ICCM, 10, Whistler, Canada, 1995, pp. I:245–252. Cui, W. and Wisnom, M.R., A combined stress-based and fracture mechanics based model for predicting delamination in composites. Composites, 1993, 24, 467–474. Robinson, P., Javidrad, F. and Hitching, D., Finite element modelling of delamination growth in the DCB and edge delaminated DCB specimens. Composite Structures, 1995, 32, 275–285. Krger, R., Knig, M. and Schneider, T., Computation of local energy release rates along straight and curved delamination fronts of laminated DCB and ENG specimens. In 34th AIAA/ASMG/ASCE/AHS/ ASC SSDM Conference, 1993, pp. 1332–1342. Mi, Y., Crisfield, M.A., Hellweg, H.-B. and Davies, G.A.O., Mixed mode delamination via interface elements. Phil. Trans. Roy. Soc., 1996, submitted. Rybicki, E.F. and Kanninen, M.F., A finite element calculation of stress intensity factors by a modified crack enclosure integral. Engineering Fracture Mechanics, 1977, 9, 931–938. Hellen, T.K., On the method of virtual crack extension. Int. J. Num. Meth. Engng, 1975, 9(1), 187–208. Schellekens, J.C.J. and de Borst, R., Numerical simulation of free edge delamination in graphite/epoxy specimen under uniaxial tension. International Conference on Composite Structures, ed. I.H. Marshall, Paisley, Scotland, 1991. Elsevier Science, Oxford. Hashemi, S., Kinloch, A.J. and Williams, J.G., The analysis of interlaminar fracture in uniaxial fibre-polymer composites. Proceedings of Royal Society, A427: pp. 173–198, London, 1990. Williams, J.G., The fracture mechanics of delamination tests. Journal of Strain Analysis, 1989, 24(4), 207–214. Davies, G.A.O. and Zhang, X., Impact damage prediction in carbon composite structures. International Journal of Engineering, 1995, 16(1), 149–170. Davies, G.A.O., Zang, X., Zhou, G. and Watson, S., Numerical modelling of impact damage. Composites, 1994, 25(5), 342–350. Robinson, R., Solutions for energy release rates of some special cases of delaminated plates. Composites, in preparation.
765
Composites Part A 28A (1997) 757–765 q 1997 Published by Elsevier Science Limited Printed in Great Britain. All rights reserved 1359-835X/97/$17.00
Shear driven delamination propagation in two dimensions
G. A. O. Davies*, P. Robinson, J. Robson and D. Eady Department of Aeronautics, Imperial College of Science, Technology and Medicine, Prince Consort Road, London SW7 2BY, UK (Received 13 January 1997; revised 20 January 1997)
This paper investigates, theoretically and experimentally, the propagation of delamination in shear mode II between internal ply surfaces in a laminated composite plate. By considering two special cases the problem is approximated as an axisymmetric model amenable to analytical solutions. In one case the circular crack front expands and in the other it contracts. It is shown clearly that only a (linear) fracture mechanics approach will explain the mode and threshold of failure, but that a strength-based criterion can be used to initiate propagation without the need for an initial flaw. Apart from explaining some curious phenomena, the simple axisymmetric analytical solutions may be used as benchmarks for testing the ability of finite element models to predict plate delamination when the crack front is curved (at present only simple beam models are available for the standard DCB and ENF tests, but even these tests do not really develop straight crack fronts). q 1997 Published by Elsevier Science Limited. (Keywords: delamination; mode II; finite element simulation)
INTRODUCTION Laminated carbon-fibre epoxy composites have been a success story whenever high specific strength and stiffness is needed. Military and civil aircraft, and racing cars are clear examples, although the total yearly production of carbon composite structures owes far more to sporting goods than transport vehicles. Nevertheless, the use of composites is still on the increase in aerospace for example, where the manufacturing costs are coming down by a combination of automated production methods and far fewer component parts compared with metal assemblies. It has been estimated that a new generation civil aircraft wing of carbon composites can now be made which is 20% cheaper, 20% lighter, and have an infinite fatigue life compared with the metal equivalent. The main weakness of laminated composites is still the very low through-thickness strength which can be less than one fortieth of the in-plane values. Furthermore both the fibres and resin in all composite systems are brittle when compared to ductile metals, so even small stress concentrations can initiate a very local failure which may then propagate in an unstable fashion. Obvious weaknesses and potential failure sources will be between bonded components or at any form of discontinuity, in load or geometry, but another hidden menace is the interior delamination due to shear stresses which will likely be a maximum in the mid-thickness region. Such * To whom correspondence should be addressed
internal delamination may not be serious for laminates loaded by in-plane tension but they are potentially disastrous for plates and shells designed to resist compression. It is for such structures, where the specific stiffness is the design criterion, that carbon composites come into their own. The history of full-scale composite structure testing has not been altogether a happy one, and a significant number of aerospace structures have failed to meet their design load due to premature through-thickness induced failure. Such weaknesses have to be overcome by either improved resin properties, through-thickness reinforcement, or more importantly by better detailed design. It is therefore vital that designers have the analytical tools available for predicting the debilitating effect of through-thickness stresses and for evaluating them in the first place. Commercial finite element codes have become very user-friendly, are equipped with better composite models, and the modern workstation is giving practitioners the capability of modelling highly complex stress fields is in 3-D regions. It is important therefore to know what to do with the answers. This paper addresses just one source of failure driven predominantly by interlaminar shear, and the commonest situation in plates where the delamination front is curved to enclose a 2-D delaminated area.
DELAMINATION FAILURE CRITERIA It is well known that using simple strength criteria, possibly couched in terms of strains, can be misleading for brittle
757
Shear driven delamination propagation in two dimensions: G. A. O. Davies et al. materials. This is particularly the case in the presence of rapidly varying stress fields. A fracture mechanics or a damage mechanics approach may therefore be a better model. However, strength criteria have traditionally been the most widely used, and the most common test for interlaminar strength is still the 3-point bending test using a short beam1 which does give reasonably consistent values for the interlaminar shear strength. One reason for its success is the fact that the rig is usually soft loading and the shear force is constant along the beam; so that once a failure initiates, the delamination is virtually instant. Even so, the initiation of failure in the complex stress field, local to the load and supports, is a challenge to model2. Any failure in a brittle material has to capture the physics properly. Stress concentrations will not be relieved by plastic flow and yet stress concentrations can be exceptionally high in laminated composites, even in the absence of clear stress raisers—for example edge effects. Moreover, strength-based criteria alone will not scale; geometrically similar structures of different size and with identical stress fields will fail at different stress thresholds. Fracture mechanics will often resolve this dilemma since it does directly address the physics of crack propagation and the balance between the structure’s strain energy release and the energy absorbed by the movement of the advancing crack front. The nature of the stress field into which that crack propagates will also need to be correctly modelled. However, fracture mechanics needs the insertion of a flaw before the prediction of crack propagation can take place and the energy release rate (or stress intensity factor) can be quite sensitive to the size of the flaw length l. Classical linear elastic fracture mechanics in isotropic materials shows strength to be proportional to l ¹1/2. If we assume that all composite structures contain natural flaws, and then use a fracture criteria to explain experimentally found strengths, we find the initial flaw to be of the order of millimetres,3 whereas microscopic examination fails to reveal flaws of such sizes4. Nevertheless, we need separately to emulate initiation and propagation (stable or unstable) consequently and the strength criteria would seem logical before the crack appears, and a fracture criterion thereafter. This strategy has been suggested by many workers5–8 and seems to work reasonably well, the main problem being the expense and the monitoring of the crack front. By far the most common approach to compare numerical models and experimental results has been to choose 1-D tests such as the double cantilever beam (mode I), the end notch flexure test (mode II) and the short beam under 3 or 4-point bending. Even these supposedly simple cases are not strictly 1-D and the crack front can develop a curved profile at the plate free edges9,10. The aim of the fracture approach is to evaluate numerically or analytically the energy release rate G when the crack front is allowed to move a small amount, and then compare this with the known critical value. The critical value may be for a pure opening mode I (G IC), for a pure shear mode II (G IIC) or for a mixed mode11. The usual way of finding G is by virtual crack closure12 or virtual crack extension13 and hope that small crack perturbations are self similar. It appears that quite small front perturbations, using
758
ultra-fine or adapting finite element meshes, can be needed, although several workers10 have successfully propagated curved and developing crack fronts along which the perturbations are chosen to ensure propagation only when G rises to its critical value. A completely different approach, which avoids tracking a crack front completely, is the use of interface elements14 in which the interface stiffness has a loading portion which rises to a ‘strength’ and decays to zero so that the area under the load–displacement curve is equal to G C. This strategy therefore models initiation and propagation and does not have a true discrete front advancing through the mesh. Unfortunately, iterative solutions of structures with softening can lead to numerical problems which need a fair amount of computing effort to overcome11. However it is not the purpose here to debate these strategies but to validate the initiation-propagation strategy for 2-D problems, where both the energy release rate and the crack front curvature are changing simultaneously.
1-D DELAMINATION IN BEAMS The most tested benchmark cases are the standard tests for finding mode I and II fracture toughness. We concentrate on the shear mode II, or Edge Notch Flexure (ENF) and End Load Split (ELS) tests15 (see later in Figure 10). In the latter case if the two arms of the split cantilever are identically balanced and the equivalent flexural modulus of the chosen laminate is E, then the applied tip load P is shared equally and produces a pure shear field at the end of the crack length a. Simple beam theory gives the tip deflection in terms of a 3 and the energy release rate, after allowing a perturbation da, as: G¼
9P2 a2 4B2 Eh2
(1)
where B is the beam width and h the depth of each split half. This energy is solely due to bending stresses (although it can be extended to include finite rotation16) and the bending stresses increase linearly away from the applied load so that the crack is advancing into an increasing stress field. The energy release therefore increases with a in equation (1) and the delamination is therefore unstable under the constant dead load P. This test is an accepted standard but the theoretical model is 1-D, and the crack front does not change in length. We now turn to two examples where this is not so.
2-D DELAMINATION IN THIN PLATES To keep the analysis simple, and free from further complications in numerical strategies, we pose two axisymmetric problems with very different internal shear fields. There is of course no such structure as a laminated composite plate of constant thickness and having axisymmetrical properties, so this is our first approximation. However, the tests described will all deploy quasi-isotropic lay-ups with a large number of plies up to 32 for thicker plates so the structural response will be close to axisymmetrical. The first case is
Shear driven delamination propagation in two dimensions: G. A. O. Davies et al. shear strength (t u) then suppose the laminate has a large number of plies in the quasi-isotropic lay-up, none of them blocked together. The maximum shear stress will then be approximately 1.5 3 the mean, and we can predict failure when: 1:5
P ¼ tu 2prt
and the circular delamination area will then be given by A ¼ pr2 ¼
Figure 1 Variations in shear and peeling stresses away from a point load
that of a single concentrated force where the far-field through-thickness stress decays like 1/r and the second case a uniform pressure, reacted by a far-field boundary, where the shear stress increases like r. Concentrated force It will be helpful to consider the nature of the mean internal stresses j zz and t rz as shown in Figure 1. The farfield shear stress t rz is simply given by P/2prt until we approach the axis of symmetry under the load where, by symmetry, it must be zero. Clearly the stress and strain field immediately under the contact force is highly complex and nonlinear, but the stresses will approach the far field values in a transition region of order t as indicated in Figure 1. From the equilibrium equations it is found that the throughthickness peel stress j zz is proportional to the derivative of the shear and is sketched in the figure with a compressive zone under the applied force and a tensile zone further away, in which numerical models have shown the stress to be small in magnitude. The zone vulnerable to delamination is therefore away from the load when the compressive stresses change to tensile and the shear is still high. Highly detailed finite element models confirm this ‘average’ description and indeed tests show sometimes a delaminated annulus surrounding an intact centre17. If we wish to use a strength criteria, based on interlaminar
2
9 P 2 t 16pt u
(2)
Because the shear stresses decrease to zero in the vicinity of the origin, we do not expect equation (2) to be valid for a distance of, say, r ¼ t. A somewhat arbitrary bound on this model is therefore noted as A ¼ pt 2. As it happens a very large number of low velocity impact tests, on an equally large variety of plate sizes and thicknesses, were available for testing equation (2). The material was HTA/6376 for all plates. These tests have been summarised in Davies and Zhang17 but for completion Figure 2 has been reproduced here to see if equation (2) is a valid model. In low velocity impact, the force generated is many times greater than the distributed plate inertia forces, so the above static model is fair. However, we admit that delamination is far from a single circular patch. The damage is not confirmed to a single plane—in fact the middle plane of a symmetrical stacking sequence is not even a true interface. Instead the delamination consist of elongated oblongs in the direction of the ply fibres at the delamination lower surface, and these delamination may occur on four or more interfaces in the central region of the plate thickness. However the envelop, as revealed in ultrasonic scans, is nearly a perfect circle17. It appears that the envelop of these oblong fronts does advance, as a circle for quasi-isotropic laminates, with increasing load. With these reservations Figure 2 shows the application of equation (2) to a series of impacted plates of many sizes and having 8, 16, and 32 plies. An interlaminar shear strength of 50 MPa for HTA/6376 composite is a lower bound to take into account the fact that delamination occurs between different ply orientations. We see that the three strength curves, even with the lower bound imposed, do not reflect the behaviour of the damage. For the thinnest plates of 1 mm the map is reasonably up to a force of 1000 N, but for the thickest plate equation (2) completely misses the very rapid development of damage which occurs at an impact force of approximately 5500 N. In fact it was not possible in these tests to control the damage near this load, the development of delamination was truly unstable, and suggested therefore a fracture mechanism. If we therefore turn to Figure 3, taken from Davies et al.18, we can evaluate energy release rate, around the boundary of a single circular delamination, in the same way as equation (1) using isotropic plate theory for the deformation. If the delamination did not exist in the plate zone bounded by A and A9B9 we can easily find the deflection of the load P due solely to the strains in the region ABB9A9 and taking the normals like AA9 and BB9 as fixed in space. Denote this
759
Shear driven delamination propagation in two dimensions: G. A. O. Davies et al.
Figure 2
Strength-based predictions for delamination area. (X denotes lower bound)
therefore the true extra deflection of P due to the insertion of the delamination, and the net work available to drive the crack is given by: 1 9Pa2 (1 ¹ v2 ) U ¼ P(du ¹ dd ) ¼ 2 64pEh3 Notice that this energy is a function of a 2, as for the 1-D beam equation (1), but the crack front length increases with a, so that when we equate: G·2pa da ¼
]U da ]a
we find G¼
Figure 3 Simplified model of internal delamination for a quasi-isotropic
If there is a fraction-free delamination of radius a then the two separate halves move together and share the load P equally, whence the ‘damaged’ deflection is this time given by: P dd ¼ 3 a2 (1 ¹ v2 )4ph3 2 The normals AA9 do of course rotate, but the rotation is due to the strains in the undamaged region outside ABB9A9, and is the same for both cases. The difference d d ¹ d u is
760
(3)
Equation (3) is remarkable in that it is independent of the size of the delamination, a. (In fact, in view of the difficulties in verifying crack front positions in tests, it might be a candidate for replacing the ELS test!). It does successfully explain the curious behaviour of the damage maps in Figure 2, if we equate equation (3) to the critical value of G IIC and invert. The critical force level
‘undamaged’ deflection as: du ¼ 3Pa2 (1 ¹ v2 )=4p(2h)3
9P2 (1 ¹ v2 ) 64p2 Eh3
P2c ¼
8p2 E(2h)3 GIIC 9(1 ¹ v2 )
(4)
and is also independent of crack size. Thus there is a genuine threshold force, below which no damage will propagate, but when the threshold is reached the amount of damage is indeterminate if the threshold force is kept constant. There is nothing to be gained by trying to use fracture analysis to evaluate the extent of damage. It also explains why there is so much scatter in damage maps. It is not due to variable non-reproducible composite fibre and resin properties, as is often suggested.
Shear driven delamination propagation in two dimensions: G. A. O. Davies et al.
Figure 4
Fracture-based predictions for force threshold
In Davies and Zhang17, equation (4) was used with an average value of G IIC of 0.8 N/mm to account for the differing interfaces at which delamination occurs. Figure 4 shows that the thresholds are predicted reasonably well. For the set of thickest two plates the strength-based initiation is not activated at all, and the lower bound from Fig. 3 is believable. The unstable fracture threshold is however predicted accurately. For the thinnest set of plates, the lower bound does not now inhibit the strength-based curve as much and the tests do show the parabolic rising curve until an unstable growth at about 1000 N. The fracture-predicted threshold value on this occasion is a little on the low side.
The point to note is that, although the shear stress decreases away from the outer boundary, so does the perimeter of an advancing crack front. It is shown19 that as the crack advances from the outer edge of the circular plate of radius R, by a radial distance l, the mode II energy release rate is GII ¼
9 p2 (1 ¹ v2 ) [(2R ¹ l)l]2 32 Et3
Inverting this equation there should be a critical pressure, given in terms of G IIC, as
4 2GIIC Et3 1 Pc ¼ 2 3 (1 ¹ v ) (2R ¹ l)l
Uniform pressure The other case of a plate loaded by uniform pressure, and needing therefore an outer edge support, has completely different stress field, again amenable to simple theory. If the applied normal pressure is p, the average shear stress at any radius distance r from the centre, will be given by: ppr2 =2prt and allowing for the maximum stress t to be 1.5 3 the mean, we have: t¼
3 pr 4 t
(5)
Equation (5) can be used to estimate the radius of the inner boundary of delamination by putting t equal to the interlaminar shear strength. As the shear increases linearly the delamination will extend from the outer support to the front but will apparently be a stable delamination since the stress decreases from the delamination front into the undamaged central region. The true test of stability needs an energy release approach. In this case the work is done by a distributed pressure, but it is amenable to an energy release analysis similar to the previous example19.
(6)
(7)
Unlike the previous example this pressure loading does depend on the delamination length l and for small l p R, it varies inversely like 1/l. As l increases the critical pressure in equation (7) decreases, indicating that the delamination is unstable, and not as the strength criteria would have us believe. This is not too dissimilar from in-plane classical fracture phenomena where the applied critical tension varies like (crack length) ¹1/2, so we have returned to the dilemma of selection an initial flaw. We therefore have conducted a series of experimental tests to investigate the true nature of the failure mode, initiation, and development of the crack fronts.
EXPERIMENTAL DETAILS To pick up possible scale effects, two sizes of circular plates were chosen, 130 m diameter and 270 mm diameter. The plates, this time made of T800/924 carbon epoxy, were held with the edges clamped and an artificial delamination using 0.0125 mm thick Teflon insert was placed in the mid-plane during fabrication. In a symmetrical lay-up the mid-plane is not a true interface, so to represent a natural delamination
761
Shear driven delamination propagation in two dimensions: G. A. O. Davies et al.
Figure 5
Figure 6
Test rig for applying uniform pressure
History of displacement and pressure in a typical test
source the adjacent plies were chosen to be 6 45 by using the following 32-ply (nominal 4 mm) lay-up [A= ¹ A=A= ¹ A] where A ¼ 45/90/ ¹ 45/0/0/ ¹ 45/90/45.
762
This lay-up is balanced with no twist–bend coupling and was perfectly flat when removed from the autoclave. Delamination widths of 5, 10, 15, 20, and 25 mm were chosen. The plates were clamped by bolts and gasket as shown in Figure 5, which sealed the container of water to be pressurised using a 200 bar hand-pump. The water pressure and the plate’s centre deflection were monitored continuously during the manual loading. The loading was increased at as near constant rate as possible until the delamination was clearly heard to propagate. ‘Failure’ was instant and unambiguous, and the pressure dropped immediately with very little change in plate deflection since water is virtually incompressible. Figure 6 shows a typical history. After failure, three post-mortem diagnoses were used. Each plate was ultrasonically C-scanned to confirm the unstable fracture, and Figure 7 shows a typical output. The original is more clearly colour coded, to reflect the time-offlight to the delamination depths, and shows that the artificial (5 mm) delamination has propagated into the plate centre over the regions subtending roughly 10 to 2 o’clock
Shear driven delamination propagation in two dimensions: G. A. O. Davies et al.
Figure 7 Final delaminated zones starting from 5 mm around the edge. The grey scale is the mid-plane
and 3 to 9 o’clock. Away from the axes of symmetry the 6 45 interface becomes 0/90 to an advancing circular front, and here the delamination switched to the weaker 90/45 interface. Nevertheless, the unstable front is almost axisymmetrical and probably would have advanced to the place centre had the pressure loading between sustainable. Fractographic examination of the failure surfaces showed characteristic cusps indicating a dominant mode II fracture mostly in the radial direction. A very small region of mode I peel was visible at the edge of the Teflon insert. Finally it was predicted that bending stresses around the clamped edges would exceed 2000 MPa for the large plates. On taking sections the polished micrographs confirmed extensive compressive microbuckling near the top surface and cracking near the lower (Figure 8). The consequent degradation of the plate’s bending stiffness would modify the clamped boundary condition, but fortunately this does not affect the simple theoretical predictions of equations (5) and (7). Turning to these equations we need to select an interlaminar shear strength, a fracture toughness, and a flexural modulus. Unfortunately, all three vary with orientation around the plate. There are also different values of fracture toughness for initiation and propagation. We therefore selected upper and lower bounds of: Interlaminar Shear Strength: 120 and 80 MPa G IIC: 1.1 and 0.55 N/mm E: 60 and 44 GPa These bounds are rather widely spaced, although as indicated in equation (7) the critical pressure depends on the square root of the product of toughness and modulus. Nevertheless the results plotted in Figure 9 do appear conclusive for both plate sizes. The shear-strength based criteria does predict the initiation of failure for the undamaged plates, but becomes progressively inaccurate as the initial delamination size increases. The fracture-based
Figure 8 Micrograph near top and bottom surfaces under high bending strains
criterion is quite respectable in spite of the simplifying nature of our theory. Thus the experiments do confirm the nature, mode, and critical pressure of the theoretical failure model. The test data clearly switches from strength to fracture even though the transition is not quite sharp as the two separate theories would indicate. In this case, rather like the 1-D short beam test, the shear strain does not vary drastically for a small distance from the plate edge, so there is no need to explain why an undamaged plate needs an initial flaw for the fracture criteria to be triggered. It is permissible to use the strength predictor until the delamination length causes the fracture criteria to override it. CONCLUSIONS • Without introducing the complexities of a numerical analysis it has been possible to select two (quasi) axisymmetrical problems of plate bending where the failure mode is delamination, driven by mode II shear, and the crack front envelope is nearly circular. • In both cases a strength criterion can be used to initiate a delamination, but the energy release model is needed to
763
Shear driven delamination propagation in two dimensions: G. A. O. Davies et al.
Figure 9
(a) Small plate failure. (b) Large plate failure
predict the subsequent growth and the stability of that growth. • A numerical simulation using interface elements should be successful if the load/deflection law increases to a maximum strength but thereafter decreases so that there is total energy release equal to G IIC. • A fracture based virtual crack closure (or extension) strategy should simulate the advancing crack front, but it is suggested that the initial node releasing may be based on interlaminar strength until (and if) the fracture analysis indicates growth at a lower load than the strength based criterion. • Finally we note that there is much current research aimed at using Finite Element models to initiate and propagate delamination. In reality most crack fronts are curved,
764
even for the simple beam tests illustrated in Figure 10, and the F.E. community has to validate the numerical models with experimental observation involving all the usual uncertainties of strength, modules, toughness, and the position of the crack front. We therefore offer these two axisymmetrical solutions as a benchmark for curved crack front propagation.
ACKNOWLEDGEMENTS The authors would like to thank the DRA (Structural Materials Centre, Farnborough) for supporting projects on low velocity and ballistic impact which revealed the need for these basic failure studies.
Shear driven delamination propagation in two dimensions: G. A. O. Davies et al. 6. 7. 8. 9. 10.
11. 12.
Figure 10 Standard ENF and ELS tests for mode II fracture toughness
13. 14.
REFERENCES 1.
2. 3. 4. 5.
Curtis, P.T., ed., Crag test methods for the measurement of the engineering properties of fibre-reinforced plastics. Royal Aerospace Establishment, Tech. Rep. 88012, Farnborough, UK, February 1988. Wisnom, M.R., Cui, W.C. and Jone, M., Failure mechanisms in 3 and 4-point short beam bending tests of u-d glass/epoxy. Journal of Strain Analysis, 1992, 27(4), 235–243. Davies, G.A.O. and Cui de, Yu., Three dimensional edge effects in composite plates. Computers and Structures, 1980, 31(1), 11–16. Hull, D., An Introduction to Composite Materials. Cambridge University Press, 1981. Bostaph, G.M. and Elber, W., A fracture mechanics analysis for delamination growth during impact on composite plates. ASME Winter Conference, 1982.
15. 16. 17. 18. 19.
Wisnom, M.R., Interface element approach to modelling mode I and II fracture. In ICCM, 10, Whistler, Canada, 1995, pp. I:77–83. Minguet, P.J. and O’Brien T.K., Analysis of composite skin/stringer bond failure using a strain energy release rate approach. In ICCM, 10, Whistler, Canada, 1995, pp. I:245–252. Cui, W. and Wisnom, M.R., A combined stress-based and fracture mechanics based model for predicting delamination in composites. Composites, 1993, 24, 467–474. Robinson, P., Javidrad, F. and Hitching, D., Finite element modelling of delamination growth in the DCB and edge delaminated DCB specimens. Composite Structures, 1995, 32, 275–285. Krger, R., Knig, M. and Schneider, T., Computation of local energy release rates along straight and curved delamination fronts of laminated DCB and ENG specimens. In 34th AIAA/ASMG/ASCE/AHS/ ASC SSDM Conference, 1993, pp. 1332–1342. Mi, Y., Crisfield, M.A., Hellweg, H.-B. and Davies, G.A.O., Mixed mode delamination via interface elements. Phil. Trans. Roy. Soc., 1996, submitted. Rybicki, E.F. and Kanninen, M.F., A finite element calculation of stress intensity factors by a modified crack enclosure integral. Engineering Fracture Mechanics, 1977, 9, 931–938. Hellen, T.K., On the method of virtual crack extension. Int. J. Num. Meth. Engng, 1975, 9(1), 187–208. Schellekens, J.C.J. and de Borst, R., Numerical simulation of free edge delamination in graphite/epoxy specimen under uniaxial tension. International Conference on Composite Structures, ed. I.H. Marshall, Paisley, Scotland, 1991. Elsevier Science, Oxford. Hashemi, S., Kinloch, A.J. and Williams, J.G., The analysis of interlaminar fracture in uniaxial fibre-polymer composites. Proceedings of Royal Society, A427: pp. 173–198, London, 1990. Williams, J.G., The fracture mechanics of delamination tests. Journal of Strain Analysis, 1989, 24(4), 207–214. Davies, G.A.O. and Zhang, X., Impact damage prediction in carbon composite structures. International Journal of Engineering, 1995, 16(1), 149–170. Davies, G.A.O., Zang, X., Zhou, G. and Watson, S., Numerical modelling of impact damage. Composites, 1994, 25(5), 342–350. Robinson, R., Solutions for energy release rates of some special cases of delaminated plates. Composites, in preparation.
765