Jtam-A4.dvi JOURNAL OF THEORETICAL AND APPLIED MECHANICS 52, 2, pp. 301-312, Warsaw 2014 NUMERICAL AND EXPERIMENTAL INVESTIGATIONS OF EMBEDDED DELAMINATION GROWTH CAUSED BY COMPRESSIVE LOADING Piotr Bajurko Institute of Aviation, Warsaw, Poland; e-mail: piotr.bajurko@ilot.edu.pl Piotr Czarnocki Warsaw University of Technology, Institute of Aeronautics and Applied Mechanics, Warsaw, Poland e-mail: pecz@meil.pw.edu.pl Thepaperdealswith growthanalysis of initially circulardelaminations embedded in carbon- epoxy laminate plates subjected to compressive loading. Three different reinforcement lay- ups yielding different elastic laminate properties were considered. The numerical results were supplemented with experimental ones. The reasonably good agreement between the numerical predictions and experimental results was found. It was shown that variation in elastic properties of sub-laminates separated by delaminations significantly affected the way the delaminations propagated. Keywords: buckling, driven, delaminations 1. Introduction For strength and stiffness reasons, composite thin-walled airframe parts consist of a large num- ber of reinforcement layers forming sub-laminates of the same relative fibre orientation. The sub-laminates, in turn, differ in the relative reinforcement orientation and such a difference be- tween the adjacent sub-laminates is equal to 45◦ inmost cases. Due tomanufacturing processes consisting in the layer-by-layer placement and their consolidation, typical manufacturing de- fects result from the air entrapped between the adjacent reinforcement layers. The entrapped air forms embedded delaminations that can grow when buckle under compressive loading. On- ce such a defect has been detected in the course of quality check, the decision must be made whether the defected airframe part must be scrapped or it can be repaired or, perhaps, left as it is. The answer is not obvious and depends on the delamination growth ability which must be investigated. It can be done with the use of the Linear Fracture Mechanics (LFM) tools taking advantage of appropriate fracture criteria by comparing certain combinations of expec- ted (calculated) values of the strain EnergyRelease Rate (SERR) components GI,GII and GIII against corresponding critical ones being thematerial constants determined experimentally. The former have been addressed in a number of papers in which various delamination models were presented. Early publications of Kacchanov (1976), Chai et al. (1981), Bolotin (1984, 1988) and Hutchinson et al. (1987) dealt with 2Dmodels which offered closed analytical formulas butwere of limited practical applicability due to simplifications in the model geometry. The models co- uld be applied for beam delaminations or throughout delaminations in plates of infinite width. First, more realistic 3D models representing embedded delaminations and focusing on critical buckling loading yielding out of plane deformation of a disbanded laminate layer were presented byShivakumarandWhitcomb (1985),Kassapoglou andHammer (1990),KimandMayer (2003). More advanced models allowing for analysis of delaminations at the post-buckling stage were presented by Konishi (1985), Chai and Babcock (1985), Chai (1990), Jane and Yin (1992), Yin and Jane (1992), Bolotin (1996). Themodel presented byKonishi (1985) allowed for estimation 302 P. Bajurko, P. Czarnocki of the stress state at certain points in the vicinity of delamination front. Chai (1985) andBolotin (1996) assumed a self-similar delamination growth (which is not realistic) and determined the distribution of total SERR values along the delamination front. Chai (1990) determined the values of GI and GII for delamination of an isotropic laminate at certain delamination front points eliminating the simplifying assumption of the self-similar delamination growth.Models of similar capability, however for laminates of orthotropic layers were developed by Jane and Yin (1992), Yin and Jane (1992) and Sheinman and Kardomateas (1997). Ritz’s method was used to determine the post-buckling delamination geometry and was combined with methods taking the advantage ofmembrane forces and bendingmoments to determine the SERRvalues. An im- portant feature and shortcoming of all Ritz’s method basedmodels consisted in the need for an a priori assumption on the buckling geometry. Reported, more advancedmodels obtained with the use of the FEMdid not offer solutions in closed forms either, however, they did not require any assumption on the buckling geometry and the way it would grow. For example, the FEmo- dels presented by Whitcomb (1990), Withcomb (1992) represented elliptical delaminations and allowed for discrete determination of all components of SERRwith the use of themodified crack closure method first proposed by Rybicki and Kanninen (1997). Whitcomb (1990) for the first time made a provision for possible initial contact between the delamination separated layers, however because of number of simplifications, the regaining of once lost contact in the course of delamination buckling geometry changes was not possible. More advanced FEmodels taking advantage of typical contact elements and allowing for variation in the contact regions occurring in the course of buckling growthwere presented byCzarnocki (2000) andRiccio et al. (2001). In early publications, the growth of fully embedded delaminations and resulting from it progres- sive laminate damage were hardly addressed due to the lack of appropriate tools to tackle the problem. Development in the Virtual Crack Closure Technique (VCCT) and cohesive element formulations provided effective tools for this purpose. The cohesive element technology is based on early ideas ofDugdale (1960) andBarenblatt (1962). Nowadays, a detailed description of this technology can be found in a number of papers or commercial FE codesmanuals e.g. ABAQUS v.6.10 Manual. The cohesive element technology is an especially suitable tool for investigation of damage caused by an interfacial failure since, in general, cohesive elements are interface ones and can be placed between the solid elements representing the structure. The failure of cohesive elements is definedwith the prescribed failure criterion, i.e. the delamination criterion resulting in separation of solid elements and can mimic the interfacial failure. Although, the use of this technology allows for separation at the prescribed plane only, nevertheless, an appropriate plane is easy for anticipation due to a laminate structure. Using this technology, a number of delami- nation tests was simulated such as delamination under Mode I, Mode II and Mixed-mode I/II loading conditions, e.g. Fan et al. (2008), Balzani and Wagner (2008) and Corona and Reedy (2011), respectively. The FE models developed for this purpose were relatively simple and re- presented pre-cracked beams or plates with through-width delaminations propagating in a near self-similar manner. Such models were of limited practical application and could serve rather for the purpose of the cohesive element technology testing. A fewmost recent papers presented more usefulmodels of fully embedded delaminations of rectangular and circular boundaries, see Wang et al. (2010), Lampani (2011) and Butler et al. (2012). Possible growth of buckling-driven fully embedded delaminations depend on the delamina- tion resistance of the interface between the delamination separated layers as well as on elastic properties of delamination created sub-laminates which, in turn, depend on the elastic proper- ties of reinforcement layers and their stocking sequence. Fibre relative orientationmismatch can affect the fracture resistance, however, for the relative reinforcement orientation of interest the effect of such a mismatch is not strong, see Pereira and deMorais (2004a,b). The research presented below was focused on the relationship between the mismatch of the elastic properties of sub-laminates caused by the delaminations and their growth under in-plane Numerical and experimental investigations of embedded delamination... 303 compressive loading. For the purpose of analysis, one assumed that the fracture toughness of the laminate under consideration did not varied with the relative fibre orientation. Such an assumption could be justified by the results presented by Pereira and de Morais (2004a,b). A plate of three different reinforcement layer stockings with initially circular delamination located in the centre of theplatewas considered,Fig. 1.Thedelaminationwas located in such away that its plane divided the laminate into two sub-laminates of different thicknesses. In each case, the thickness ratio h/H =1/6, where h and H are the thicknesses of the thinner and thicker sub- laminates, respectively. Delamination growth was modelled with the help of cohesive elements for three different stocking sequences of the reinforcement layers: [0◦4//0 ◦ 24], [0 ◦ 4//90 ◦ 20/0 ◦ 4] and [45◦/−45◦/−45◦/45◦//90◦20/45 ◦/−45◦/−45◦/45◦] where the symbol // denoted location of the delamination. Open literature on growth of such delaminations, especially for the last mentioned configu- ration, is limited and one could expected that it would be difficult to verify the obtained results against those found in the literature. For this reason, an experimental verification of the nume- rical analysis results for this configuration was carried out. For this purpose, a layered plate of [45◦/−45◦/−45◦/45◦//90◦20/45 ◦/−45◦/−45◦/45◦] reinforcement arrangement and the same geometry as in the numericalmodel wasmanufactured.An artificial delaminationwas produced by aTeflon film circular insert located between the respective reinforcement layers in the course of the lamination process. The plate was loaded under the displacement controlled conditions in a fixture recommended by ASTMD7137/D7137M to prevent global plate buckling. A set of six strain gauges located in the selected plate regions was used formonitoring the delamination buckling and delamination growth. Since the results of the preliminary FE analysis indicated that the changes in εy were amore sensitive indicator of the buckling and delamination growth than the changes in εx, all gauges were arranged in the y direction. The same analysis was used in selection of plate regions for which εy fields were almost uniform and of very low gradients, and were not affected by the delamination changes in the course of loading, e.g. location V. Strain in this location was related to the displacement of the loaded plate edge for both the actual specimen and its numerical model. Subsequently, this relation was used for comparison of the numerical results to the experimental ones. Since the method involving the use of strain gauges provided a discrete type of informa- tion limited to the points of gauge location, a C-scan and computed tomography (CT) post mortem inspections of the defected plate region were carried out to obtain more comprehensive information about the final delamination geometry. 2. Numerical analysis 2.1. Geometry of the modelled structure A rectangular carbon-epoxy laminate plate 150mm×100mm×4mm containing an embed- ded, circular 40mmdiameter delamination located in the centre of the plate below the 4th rein- forcement layer, (approximately 0.6mmbelow theplate surface),wasmodelled, Fig. 1.Theplate geometrywas chosen ensuring that for the experimental verification of the numerical results, the standard Compression after Impact (CAI) fixture could be used (ASTMD7137/D7137M). The plate consisted of 28 reinforcement layers. Each impregnated reinforcement layer was approxi- mately 0.144mm thick. Elastic properties of the single cured layer are presented in Table 1. 2.2. FE representation The geometry of the FE model, i.e. the mesh density and types and shapes of elements modelling the composite structure as well as the boundary conditions were the same for all 304 P. Bajurko, P. Czarnocki (a) (b) Strain gage Co-ordinates No. x [mm] y [mm] Front side I 0 0 II −20 0 III −35 0 IV 0 −20 V 0 −50 Back side Ia 0 0 IVa 0 −20 Fig. 1. A view to the analyzed structure-front side (a) and localization of strain gauges (b). All dimensions given in mm Table 1.Elastic properties of a single laminate layer E11 E22 =E33 G12 =G13 G23 ν12 ν23 = ν13 [MPa] [MPa] [MPa] [MPa] 128290 8760 4270 3000 0.288 0.320 three laminate structures. Themeshwas homogeneous and composed of square cuboid elements (Type117 in theMSCMARCnomenclature) of 1mm×1mmbaseFig. 2.Toavoid anydependen- cy between the delamination geometry changes andmesh geometry, a rectangular homogenous mesh was applied. The boundary conditions applied to the model reflected the way the speci- menwas supportedwhen tested with the use of CAI standard rig. The initial delamination was defined by location of nodes whoose x and y coordinates met criterion √ xx+y2 < d 2 (2.1) where d was the nominal value of the delamination diameter and x and y co-ordinates were expressed in the global coordinate system with the origin in the centre of the plate. Fig. 2. Finite element representation of the analyzed structure (a) and FE definition of the initial delamination geometry (b) This procedure resulted in replacing the circular delamination of the assumed diameter d with a polygonwith not smooth, saw-like boundary following the shape of theFEmesh, Fig. 2b. The plate was compressed in the y direction while the fixture prevented displacements in two Numerical and experimental investigations of embedded delamination... 305 other directions. In the expected delamination plane, the interface elements of type 188 with zero thickness were located. These elements were associated with a cohesive zone material of bi-linear traction-displacement relationship, Fig. 3. The in-plane size of each interface element was equal to the size of the adjacent solid elements. Thebi-linear relation between the traction t and separation δ was defined with three parameters: critical value of the SERR Gc, maximum separation δm, and separation δc, corresponding to themaximumtraction tmax, Fig. 3.A failure of the cohesive element occurred if δ > δm triggering delamination growth. In theMSCMARC formulation, for most general loading cases three relative displacement components could be considered: out-of-plane, opening displacement δn and two in-plane sliding displacements δs and δt. These displacements could be used to define the effective displacement δ= √ δ2n+β 2 1(δ 2 s + δ 2 t ) (2.2) where β1 = ts,t/tn is a coefficient used to allow for the difference in the sliding and opening displacements. To definedifferent values of Gc in shear than in tension, the shear/normal energy ratio β2 could be used. In a general state of deformation, when β1 6=1, the curve defining the effective traction versus the effective opening displacement was defined as a linear combination of the response in pure tension and pure shear. Cohesive material parameters GIc =300 ∗N/m β2 =3 tnmax =30 ∗MPa β1 =2 δIc =0.00033 ∗∗mm δIm =0.018 ∗∗mm ∗) TEBUKprogress rapport, Institute of Aviation (2011) ∗∗) calculate based on Harper (2008, 2010) Fig. 3. Constitutive law of the cohesivematerial model The FE analysis consisted of the following essential steps: 1) linear buckling analysis to determine the critical buckling load of the thinner sub-laminate in the region of delamination 2) application of small internal pressure acting on the delamination faces and simultaneous application of compressive loading exceeding by 10% the critical buckling load determi- ned by the linear buckling analysis to produce a slight out-of-plane displacement of the delamination separated layer 3) removal of the internal pressure and increasing of the compressive loading up to the pre- scribed maximum value to simulate delamination growth. 2.3. Results of FE modelling Figures 4 to 6 present contour lines for progressing delaminations embedded in [0◦4//0 ◦ 20/0 ◦ 4], [0 ◦ 4//90 ◦ 20/0 ◦ 4] and [45 ◦/−45◦/−45◦/45◦//90◦20/45 ◦/−45◦/−45◦/45◦] lami- nates respectively. Each set of the contour lines is supplemented with the corresponding load- displacement graph. Inaddition, the correspondingcontour lines andpoints of load-displacement graphs are marked with the same digits. Graphs in Fig. 7 present the compressive force versus delamination area. 306 P. Bajurko, P. Czarnocki Fig. 4. Delamination growth for 0◦4//0 ◦ 24 reinforcement arrangement: (a) delamination contour lines, (b) corresponding points of the force-displacement curve Fig. 5. Delamination growth for 0◦4//90 ◦ 20/0 ◦ 4 reinforcement arrangement: (a) delamination contour lines, (b) corresponding points of the force-displacement curve Fig. 6. Delamination growth for 45◦/−45◦/−45◦/45◦//90◦20/45 ◦/−45◦/−45◦/45◦ reinforcement arrangement: (a) delamination contour lines, (b) corresponding points of the force-displacement curve. Absolute force and displacement values are shown Fig. 7. Force versus delamination area. Open symbols represent the buckling force corresponding to the initial circular delamination Numerical and experimental investigations of embedded delamination... 307 Graphs in Fig. 8 represent the calculated and experimental force versus the longitudinal strain εy relationships. The strain values correspond to that in location V. Fig. 8. Calculated and experimental load-strain curves. Absolute force and strainsvalues are shown 3. Experimental work To gain information about the quality of the numericalmodel developed, the experiment similar to the CAI test was carried out following the recommendations of ASTM D 7137/D7137M- 05 standard. The changes in the area and shape of delamination resulting from the growth of initially circular delaminationwere investigated with the help of theCTand ultrasoundmethod providing a C-scan of the damaged region. In addition, to compare the measured strain values with the numerical predictions, a local measurement of strains in the course of loading was carried out with the use of five strain gauges in locations I, II, III, IV andV at the front side of the specimen and two strain gauges in locations Ia and IVa at the plate back side, Fig. 1a. The nominal x and y co-ordinates of the strain gauges centre locations are given in Fig. 1b. 3.1. Experimental set-up The compressive test was run under displacement control at 0.5mm/min cross head speed with the use of MTS320 testing machine equipped with a 200kN load cell calibrated for 50kN. The plate with the attached strain gauges of 3.2mm effective gauge length located as shown in Fig. 1a was secured with the help of the rig designed according to the ASTM D7137/D7137M recommendations to prevent global plate buckling. The strain changes were recordedwith 10Hz rate. The post mortem CT inspection was made with the use of the phoenix v|tome|x L 240 high-resolution microfocus computed tomography system. The C- scan was obtained with the use of the phase array method and the Olympus OmniScanMX system. 3.2. Results of experimental work The results of post mortem investigation of the final delamination geometry performedwith ultrasounds and radiography are shown in Fig. 9. The scans correspond to the fourth delami- nation contour shown in Fig. 6. For the visualization purpose, this contour and the contour of initial delamination were superimposed on C-scan and CT scan pictures. Diagrams in Figs. 10-12 represent relationships between the strains εy at locations I, II, III, IV, Ia and IVa versus the strains εy at location V. The continuous lines represent these relationships obtained from the strain gauge readings and the dashed lines represent the same relationships obtained from the FE calculations. 308 P. Bajurko, P. Czarnocki Fig. 9. C-scan pictures (a), CT scan (b). Both scans correspond to delamination contour No. 4 in Fig. 6. Dashed and continuous lines indicate the initial delamination contour and the final one determined by the use of FEmodelling, respectively Fig. 10. Strain εy at location I (a) and location II (b) versus strain εy at location V Fig. 11. Strain εy at location III (a) and location IV (b) versus strain εy at location V 4. Discussion It can be clearly seen from the comparison of diagrams in Figs. 4-6 that the growth of each delamination was affected by the elastic properties of the sub-laminates surrounding the dela- minations. For each case of the reinforcement lay-up, the delamination growth was different. A symmetrical growth relative to the xz and yz symmetry planes of the plate contour oc- curred for the first laminate only. In the case of the second laminate, unsymmetrical growth was observed with respect to the xz symmetry plane. Such growth could be attributed to the Numerical and experimental investigations of embedded delamination... 309 Fig. 12. Strain εy at location I a (a) and location IVa (b) versus strain εy at location V anti-symmetric delamination buckling mode. Additional linear buckling analysis indicated that there was only a small difference between the loads corresponding to the first symmetrical and the second anti-symmetric bucklingmodes. Although the first bucklingmode was symmetrical, the increasing loading produced a deformation change and transition of out-of-plane deforma- tion to amore stable anti-symmetric mode in which a half of delamination tended to penetrate the bottom sub-laminate and then was blocked, see the bottom delamination region in Fig. 4. Similar results were obtained by Butler et al. (2012). In such a deformation state, a section of the delamination front was free of peeling stress (Mode I loading) and, therefore, less prone to delamination growth which was reflected by a relatively small growth of the bottom delami- nation part, Fig. 4. In the case of laminate #3, the effect of shear-bending coupling produced by +45◦/−45◦/−45◦/+45◦ reinforcement lay-up could be easily noticed. The growth of this delamination was neither symmetrical relative to the xz nor yz symmetry plane. Comparison of the numerically predicted delamination contour against that detected by the ultrasonic and x-ray inspections showed that the numerical model underestimated the delamination size. It could be attributed to a) violation of the displacement controlled loading condition due to the finite stiffness of themachine kinematics loading chain and b) disregarding of the residual stress produced in the course of curing and arising from themismatch of the thermal expansion coef- ficients of the adjacent laminate layers differing by the reinforcement orientation. In the case of the former, the elastic deformation of the loading chain could result in accumulation of elastic energy during loading and its release in the course of delamination propagation, and this way in supplying an additional energy for the formation of a new interlaminar crack surface. In both the cases, additional tests would be needed to clarify the matter. There are two other factors that could cause the observed difference between the calculated and experimental delamination contours. Firstly, is was assumed that Gc values for delamina- tion initiation and propagation did not depend on the relative orientation of fibres. It could be justified because it was shown (Pereira, 2004a,b), that this effect was not strong. Also, it was shown by the same authors that values of Gc for the parallel relative fibre orientation were the lowest. Therefore, the calculated delamination contours should be larger than the experimental ones, which was not the case. Secondly, the propagation law assumed was the one offered by the FE code used and could not be the most suitable. However, one should bear in mind that the purpose of the research presented was to investigate the effect of elastic properties of dela- mination separated sub-laminates on the delamination propagation process and to isolate the effects of the mentioned factors. Values of Gcs had to be the same and, therefore, it was not vital which propagation law was applied. Graphs in Fig. 7 indicate that for the first and second laminate reinforcement lay-ups, the increase in load produced the change in sign of the dF/dA derivative from positive to negative, indicating unstable delamination growth. 310 P. Bajurko, P. Czarnocki The abrupt changes from compressive strain to tensile one in the centre of the disbonded sub-laminate (location I) and at location II (Figs. 10a and10b), respectively, could be attributed to sudden out-of-plane deformations, i.e. out-of-plane bending of the sub-laminate in the region of delamination. Due to such a deformation, the initially compressed flat external surface of the sub-laminate was extended. Analysis of the graphs, (Fig. 8 and Fig. 10a) shows that the value of εy in location V corresponding to the calculated buckling load (marked in Fig. 8) was the same as the value of εy (marked in Fig. 10a) corresponding to the abrupt change of strain in location I. Therefore, this relation allows one to conclude that this abrupt change of εy in location I (Fig. 10a) could be attributed to buckling. Further inspection of graphs in Fig. 10a indicate that the experimental value of strain εy at location V (i.e. determinedwith the use of strain gauge) corresponding to thementioned abrupt strain sign changewas higher than that determinedwith the use of FEM,whichwas unexpected. Usually, due to lack of flatness of a real plate resulting from a not perfect manufacturing pro- cess, the experimental buckling load is lower than the calculated one. The observed abnormality could be attributed to a stabilising effect of external atmospheric pressure. Initially, disbonded sub-laminate was airtight and the air pressure was not equal at its internal and external sur- faces. It resulted in tightening the thinner sub-laminate against the thicker one with the force resulting from the air pressure difference and produced an increase in the buckling load. Once the sub-laminate buckled, intralaminar microcracks were formed and the air got in cancelling sub-laminate stabilization. Progressing out-of-plane deformation indicated by the strain gauges at locations II, III and IV (Figs. 10b and 11a,b), respectively, occurred for the lower load than that predicted by the numerical analysis, as expected. Graphs in Figs. 12a and 12b present stra- ins at locations Ia and IVa resulting from overall axial compression. A good agreement between the calculated andmeasured values is clearly visible. The effect of FE mesh density on the delamination growth simulation was not addressed in the research presented above. Nevertheless, it is an important issue and is currently being investigated. 5. Conclusions The growth analysis under compression in the y direction of initially circular delamina- tions embedded in carbon-epoxy laminate plates of three different reinforcement lay-ups was carried out with the use of FEM. The following three different lay-ups yielding three dif- ferent elastic properties of laminates were considered: [0◦4//0 ◦ 20/0 ◦ 4], [0 ◦ 4//90 ◦ 20/0 ◦ 4] and [45◦/− 45◦/− 45◦/45◦//90◦20/45 ◦/− 45◦/− 45◦/45◦]. The analysis was focused on relations between the following three factors: compressive load, delamination area and contour, and ela- stic properties of sub-laminates between which the delaminations were located. It was found that the elastic properties of sub-laminates affected both the area of delaminations and their contour shapes. To produce a certain increase in the delamination area, the highest load was needed for the first lay-up and the lowest one for the third one. A symmetrical reinforcement relative to the yz plate symmetry plane resulted in symmetrical delamination growth relative to this plane, however the symmetry of reinforcement relative to the xz plane did not resulted in symmetrical delamination growth relative to this plane in the case of [0◦4//90 ◦ 20/0 ◦ 4] lay-up. For [45◦/−45◦/−45◦/45◦//90◦20/45 ◦/−45◦/−45◦/45◦] lay-up, antisymmetric growth relative to the xz and yz planes occurred. Also, it was found that depending on the lay-up, the delami- nations grew in a stable ([45◦/−45◦/−45◦/45◦//90◦20/45 ◦/−45◦/−45◦/45◦] lay-up) or unstable manner ([0◦4//0 ◦ 20/0 ◦ 4] and [0 ◦ 4//90 ◦ 20/0 ◦ 4] lay-ups) in the range of the load applied. Also, one could conclude that in the case of quality control of actual composite structures, attention should by paid not only to the size of detected delamination but also to its location, since the delamination growth ability is affected by the mechanical properties of delamination separated sublaminates. Numerical and experimental investigations of embedded delamination... 311 References 1. Balzani C., Wagner W., 2008, An interface element for the simulation of delamination in unidirectional fibre-reinforcedcomposite laminates,Engineering FractureMechanics,75, 2597-2615 2. BarenblattG., 1962,Themathematical theoryof equilibriumcracks inbrittle fracture,Advances in Applied Mechanics, 7, 55-129 3. Bolotin V.V., 1984, Delamination flaws in the structures of laminated composites,Mechanics of Composite Materials, 2, 239-255 4. Bolotin V.V., 1988, Interlaminar fracture of composites under combined loading, Mechanics of Composite Materials, 3, 410-418 5. Bolotin V.V., 1996, Delamination in composite structures: its origin, buckling growth and sta- bility,Composites: Part B, 27B, 129-145 6. Butler R., Rhead A.T., Liu W., Kontis N., 2012, Compressive strength of delaminated aerospace composites,Philosophical Transactions of the Royal Society A – Mathematical Physical and Engineering Sciences, 370, 1759-1779 7. Chai H., 1990, Three-dimensional fracture analysis of thin-film debonding, International Journal of Fracture, 46, 237-256 8. Chai H., BabcockC.D., 1985,Two-dimensionalmodelling of compressive failure in delaminated laminates, Journal of Composite Materials, 19, 67-98 9. Chai H., Babcocke C.D., Knauss W.G., 1981, One dimensional modelling of failure in lami- nated plates by delamination buckling, International Journal Solids Structures, 17, 11, 1069-1083 10. Corona E., Reedy E.D, Jr., 2011, Calculations of buckle-driven delamination using cohesive elements, International Journal of Fracture, 170, 2, 191-198 11. CzarnockiP., 2000,Effect of reinforcementarrangementondistributionofGI,GII andGIII along fronts of circular delaminations in orthotropic composite plates,Fracture of Polymers, Composites and Adhesives, 27, 49-60 12. DugdaleD.S., 1960,Yieldingof steel sheets containing slits,Journal of theMechanics andPhysics of Solids, 8, 100-104 13. Fan C., Jar P.Y.B., Cheng J.J.R., 2008, Cohesive zone with continuum damage properties for simulation of delamination development in fibre composites and failure of adhesive joints,Engine- ering Fracture Mechanics, 75, 13, 3866-3880 14. Harper P.W., Hallett S.R., 2008, Cohesive zone length in numerical simulations of composite delamination,Engineering Fracture Mechanics, 75, 4774-4792 15. Harper P.W., Hallett S.R., 2010, A fatigue degradation law for cohesive interface elements – development and application to composite materials, International Journal of Fatigue, 32, 1774- 1787 16. Hutchinson J.W., Mear M.E., Rice J.R., 1987, Crack paralleling an interface between dissi- milar materials,Transaction of the ASME, Journal of Applied Mechanics, 54, 828-832 17. Jane K.C., Yin W.-L., 1992, Refined buckling and postbuckling analysis of two-dimensional delaminations – I. Analysis and Validation, International Journal of Solids and Structures, 29, 5, 591-610 18. Kacchanov L.M., 1976, Failure of composite materials due to delamination,Mechanics of Poly- mers, 5, 918-992 19. Kassapoglou C., Hammer J., 1990, Design and analysis of composite structures withmanufac- turing flaws, Journal of the American Helicopter Society, 0ctober, 46-52 20. KimB.W.,MayerA.H., 2003, Influence of fiber direction andmixed-mode ratio ondelamination fracture toughness of carbon/epoxy laminates,Composite Science and Technology, 63, 695-713 312 P. Bajurko, P. Czarnocki 21. Konishi D.Y., 1985, A rational approach to the analysis of delaminated composite plates, [In:] Composite Structures, 3, J.H.Marshall (Edit.), 383-401 22. Lampani L., 2011, Finite element analysis of delamination of a composite component with the cohesive zone model technique, Engineering Computations: International Journal for Computer- Aided Engineering and Software, 28, 1, 30-46 23. Pereira A.B., de Morais A.B., 2004a,Mode I interlaminar fracture of carbon/epoxymultidi- rectional laminates,Composites Science and Technology, 64, 2261-2270 24. Pereira A.B., de Morais A.B., 2004b, Mode II interlaminar fracture of glass/epoxy multidi- rectional laminates,Composites Part A, 35, 265-272 25. Riccio A., ScaramuzzinoF., Perugini P., 2001,Embedded delamination growth in composite panels under compressive load,Composites Part B, 32, 209-218 26. Rybicki E.F., Kanninen M.F., 19977, A finite element calculation of stress intensity factors by amodified crack closure integral,Engineering Fracture Mechanics, 9, 931-938 27. Sheinman I., Kardomateas G.A., 1997, Energy release rate and stress intensity factors for delaminated composite laminates, International Journal of Solids and Structures, 34, 4, 451-459 28. Shivakumar K.N., Whitcomb J.D., 1985, Buckling of a sublaminate in a quasi-isotropic com- posite laminate, Journal of Composite Materials, 19, 2-18 29. WangR.G., ZhangL., ZhangJ., liuW.B.,HeX.D., 2010,Numerical analysis ofdelamination buckling andgrowth in slender laminated compositeusing cohesive elementmethod,Computational Materials Science, 50, 20-31 30. WhitcombJ.D., 1990,Mechanicsof instability relateddelaminationgrowth,CompositeMaterials. Testing and Design. ASTM STP, 1059, 215-230 31. Whitcomb J.D., 1992, Analysis of a laminate with a postbuckled embedded delamionation, inc- luding contact effects, Journal of Composite Materials, 26, 10, 1523-1535 32. Yin W.-L., Jane K.C., 1992, Refined buckling and postbuckling analysis of two-dimensional delaminations – II. Results for Anisotropic Laminates and Conclusion, International Journal of Solids and Structures, 29, 5, 611-639 Manuscript received May 31, 2013; accepted for print August 22, 2013