Strojniški vestnik - Journal of Mechanical Engineering 61(2015)9, 507-516 © 2015 Journal of Mechanical Engineering. All rights reserved. D0l:10.5545/sv-jme.2015.2521 Original Scientific Paper Received for review: 2015-02-22 Received revised form: 2015-06-02 Accepted for publication: 2015-06-16 Cohesive Zone Parameters Selection for Mode-I Prediction of Interfacial Delamination Mohsen Moslemi1* - Mohammadreza Khoshravan2 1 Young Researchers and Elite Club, Tabriz Branch Islamic Azad University, Tabriz, Iran 2 University of Tabriz, Department of Mechanical Engineering, 5166614766, Tabriz, Iran In order to determine the normal cohesive strength of composite laminates, a new test methodology was proposed. The values of cohesive zone parameters (the cohesive strength and the separation energy) for mode I interlamiar fracture of E-glass/epoxy woven fabrication were computed from a series of experimental tests. Cohesive zone model simulation based on interface finite elements was conducted. A modified form of the Park-Paulino-Roesler (PPR) traction-separation law together with a bilinear mixed-mode damage model was used to simulate the damage processes, using Abaqus cohesive elements. The numerical results were compared with experimental tests and confirmed the adequacy of normal cohesive strength. To ensure the sufficient dissipation of energy that successfully predicts delamination onset and propagation, cohesive zone length and minimum number of cohesive elements at cohesive zone length were determined. Interfacial penalty stiffness and the resistance curve of the composite specimen were also computed. The results show that the modified PPR model accurately simulates the fracture process zone ahead of the crack tip as compared to the bilinear model. Keywords: cohesive zone model, delamination, normal cohesive strength, finite element prediction Highlights • Interlaminar Cohesive strength of composite material was determined. • Cohesive zone modelling of delamination. • Bilinear and modified PPR damage models. • Experimental estimation of resistance curve. • Cohesive zone length estimation. 0 INTRODUCTION Delamination is a defect type that frequently occurs in laminated composite materials, described as the separation of a layer or group of layers from their adjacent ones, due to out-of-plane shear loads. Delamination usually initiates from discontinuities, such as matrix cracks and free edges, or manufacturing faults, such as embedded defects and machining process such as water-jet cutting [1]. Therefore, it is important to detect and analyse the progressive growth of delamination in order to predict the performance and to improve reliable and safe designs. Ultrasonic C-scan is one of the most efficient technics for tracing delamination in laminated composite materials [2]. Mode-I interface cracking is one of the most frequently modes of delamination in composite layered materials that is due to the loss of cohesion between layers of material. The delamination resistance of polymer matrix composites has been extensively investigated in the framework of fracture mechanics [3]. In order to model crack propagation, it was assumed that the delamination propagates when the available strain energy release rate is greater than or equal to a critical value, which is considered to be a mechanical parameter of the interface. Techniques such as the virtual crack closure (VCC), the J-integral, and the virtual crack extension are some of the most prevalent procedures that are used to predict delamination growth. These techniques are used to compute the distribution and components of the strain energy release rate. However, there are some difficulties when these techniques are performed using finite element codes. Another method to the numerical modelling of the delamination growth can be developed within the framework of damage mechanics. Models formulated according to damage mechanics are based on the concept of the cohesive crack model, which considers a zone of vanishing thickness ahead of the crack front in order to describe more realistically the fracture nature without the use of stress singularity. The cohesive zone model was first developed by Dugdale [4], who demonstrated the concept that stresses in the material are confined by the yield stress and that a thin plastic zone is generated in front of the crack tip. Barenblatt [5] proposed a cohesive zone concept to study the fracture nature of brittle materials and to introduce a separation mechanism at the atomic scale in order to describe the real separation of materials, and to remove stress singularity at the crack tip. Hillerborg et al. [6] suggested a model similar to the *Corr. Author's Address: Young Researchers And Elite Club, Tabriz Branch Islamic Azad University, Tabriz, Iran, m.moslemi@tabrizu.ac.ir 507 Strojniski vestnik - Journal of Mechanical Engineering 61(2015)9, 507-516 Barenblatt model. Their model allowed for existing cracks to grow as well as the initiation of new cracks. A cohesive zone model was frequently used to the model fracture analysis of a different variety of materials such as metals, polymers, concrete [7], ceramics, functionally graded materials [8], and fibre-reinforced materials [9]; its range of applications will continue to expand. An important issue in conjunction with the use of the cohesive zone model is the specification of the traction-separation law. In particular, the related fracture parameters, such as cohesive strength and fracture toughness, as well as the shape of the traction-separation law, must be determined. In the case of the traction-separation law, there are some models in the literature that can be used. For example, traction-separation based on an exponential form, a trapezoidal form and the bilinear form of Zhang and Paulino [8] are traction-separation laws that have been widely adopted. Since there are no available standardized tests and due to the existence of some difficulties corresponding to the direct measurement of the theses parameters, in most cases they are determined by the comparison of a measured fracture property with numerical predictions based on an idealized cohesive zone model, cf. [7], [10] and [11]. However, cohesive strength and fracture toughness are found to have higher importance in comparison to the specific traction-separation shape that chosen for the cohesive zone modelling. Turon et al. [11] proposed a methodology to estimate the constitutive parameters for the finite element simulation of progressive delamination using a bilinear cohesive zone model. According to their methodology, the cohesive strength is proportional to the length of the cohesive zone, and this parameter should be modified with respect to the cohesive zone length. It suggests that as the cohesive strength decreases, the length of the cohesive zone should be artificially lengthened to ensure that it spans enough elements of a given size. Yuan and Li [12] investigated the effects of the cohesive law on ductile crack propagation and recommended obtaining realistic computational results; the cohesive law must be defined with proper parameters. The objective of this research is to provide a suitable methodology for the fracture characterization of delamination under pure mode I loading. In this paper, a simple structure test is proposed to compute the normal interlaminar cohesive strength of composite laminates. The values for the critical mode I strain energy release rate (GIc) and mode-I cohesive strength (oIc) were computed at room temperature. These parameters and assumed damage models, including modified PPR model [13] and triangular traction law, were used with the Abaqus COH3D8 cohesive element to simulate bond line failure in structures made from E-glass/Epoxy specimen. The resistance curve of the composite specimen computed using the experimental method and a guideline methodology was proposed for selection of mode I cohesive zone length and the minimum required number of element in the cohesive zone length to obtain successful prediction of the delamination onset and propagation. 1 COHESIVE ZONE MODEL THEORY Cohesive damage zone models relate traction to separation at an interface where a crack may initiate. Crack initiation is related to the cohesive strength, i.e., the maximum traction on the traction-separation law. When the area under the traction-separation law reaches the fracture toughness, the traction declines to zero and new crack surfaces are generated. This phenomenon is shown in Fig. 1 for two different types of traction-separation law. Damage a) onset Pure mode model A ""TRxG,, M ¡g/x -> i— |\/Mixed-modc /\model -o-T+A <5m , states that the compressive stress does not contribute to crack onset. Whenever the damage initiation criterion of Eq. (3) for the cohesive element is satisfied, the stiffness of the element is declined according to the corresponding constitutive law that is illustrated in Fig. 1. 1.2 Damage Propagation Criterion After damage initiation, the softening procedure occurs. This procedure is governed by the corresponding traction separation law as below. 1.2.1 Bilinear Traction-Separation Model As shown in Fig. 1a, in the mixed-mode constitutive law, omc and c>mc represent the cohesive strength and critical separation (separation at which crack initiates), respectively. Moreover, ¿mu refers to the separation at which failure occurs. The softening relation of cohesive elements that are governed by bilinear constitutive law can be expressed as: ct, = (1 - d) Ki8i, (4) where d is a variable that relates damage condition, which has the magnitude d = 0 when the interface -I Prediction of Interfacial Delamination 509 Strojniski vestnik - Journal of Mechanical Engineering 61(2015)9, 507-516 is undamaged, and the magnitude d = 1 when the interface is completely fractured. In numerical analysis of damage evolution, there are two crack evolution criteria: energy- and displacement-based. In the analysis presented herein, the energy-based Benzeggagh and Kenane (BK) [14] damage evolution criterion was used, as expressed in Eq. (5). G = Glc + (GIIC - GIC) G G (5) In Eq. (5), n refereed to the BK material parameter, GShear = GII +GIII and GT = GI + GII + GIII- 1.2.2 Modified PPR Model In this model, the traction force in the interface is obtained by differentiating a potential function with respect to the interface separation. Fig. 1b shows the typical modified PPR traction-separation law of cohesive zone modelling. Park et al. proposed this potential based model for mixed mode fracture [15]. This law can be used for a wide variety of ductile and brittle fractures. Since it behaves nonlinearly prior to damage, it is required to develop a user defined element in ABAQUS for this model. Bhattacharjee et al. [13] developed a modified version of the PPR model for analysis of tearing in thin soft materials. In their model, as with the triangular constitutive law, a linear elastic response was assumed before damage initiation (¿< ¿c). Therefore, this modified version of the PPR model allows for a straightforward implementation of the model in the commercial finite element code ABAQUS, using tabular capability. After damage initiation, the softening relationship in this model can be expressed as [13]: 5 , w 5 w-i G = A[w(1 -— )a (- + — )w-1 - -a(1 -f ra + A )w ], o„ a 8, (6) where ô is normal displacement, and A, w, a, ôu are PPR parameters that can be determined by satisfying all boundary conditions, including: a = a at 8=8, a = 0 at 8 = 8,, ltm ^ = 0, 8^80 08 G = J ad8, (7) where G is the total strain energy release rate. By applying the above-mentioned boundary conditions, it can be expressed as: w = - a(a - 1)X2 1 -aX2 A = fC ' C = [w(1 -X)a(W + X)w-1 - a (1 - X)a-1 (w + X)w ], (8) a where X = 8 A,- (9) ¿u can be determined using the following equation as: G = - KÔ2c + f adö. n c J (10) In numerical damage simulation, the corresponding the damage variables can be defined using Eq. (4) and the corresponding traction in Eq. (6). After generating the damage variable at all the displacements, the modified PPR model can be directly implemented in ABAQUS using the tabular capability as a function of relative displacement. 1.3 Cohesive Zone Length As shown in Fig. (2), the length of the cohesive zone lcz is introduced as the distance between the crack tip and the point where the maximum cohesive traction is achieved. Fig. 2. Cohesive zone length of DCB specimen The length of the cohesive zone at the initiation of crack growth is independent of applied load and crack length. This means that the cohesive zone length at the crack initiation is a material property. There are several models in the literature that estimated the length of the cohesive zone [4] to [6], [16] and [17]. All of the models for estimation of the cohesive zone length have similar forms as mentioned in Eq. (11). lcz = RE^-, (11) where GIc is the critical strain energy release rate, oIc is the normal cohesive strength, and R is a parameter n 0 510 Moslemi, M. - Khoshravan, M. Strojniski vestnik - Journal of Mechanical Engineering 61(2015)9, 507-516 that depends on each cohesive zone model. For instance, Irwin's model [17] carried out with R = 0.31 or in Dugdale [4] and Barenblatt [5] models, the value for R is considered to be 0.4. 2 FINITE ELEMENT SIMULATION Three-dimensional finite element models are developed in ABAQUS 6.12 [18] to model the delamination onset and growth, using two different constitutive laws. The DCB model was composed of two layers of eight-nodded solid elements, C3D8R (adherent layers), which were connected with a layer of zero thickness eight-node cohesive element, COH3D8, (cohesive layer). Adherent layers were connected to the cohesive layer by surface-to-surface tie constraints. Fig. 3 illustrates the deformed shape of the DCB specimen during crack propagation. A more refined mesh near the crack tip, the outer surface of specimen and in the damage propagation region was used. Boundary conditions included applying a vertical displacement and horizontally restraining the upper and lower edge node of the arms (Fig. 3). In order to predict an accurate propagation of delamination, it is necessary to have an adequate fine finite element discretization in the cohesive zone length to achieve a good estimation of the interlaminar stress fields. When the cohesive zone is discretized by too few numbers of elements, the distribution of tractions ahead of the crack tip is not represented accurately. Thus, in order to achieve successful FEM results, it is necessary to have a minimum number of elements in the cohesive zone length. The number of cohesive elements in the cohesive zone is: N = %, (12) where le is the length of the cohesive elements in the direction of crack propagation and lcz is the cohesive zone length. There are several guidelines for the minimum number of elements in the cohesive zone length. However, this number is not well established. Fig. 3. Deformed shape of simulated DCB with 3D elements 3 EXPERIMENTAL APPROACHES 3.1 Mode I Fracture Toughness Measurement 3.1.1 The DCB Test Fig. 4a shows the geometry and dimensions of the DCB specimens. In the process of specimen fabrication, the E-glass/epoxy composite plate with a thickness of 2h = 4.2 mm was first prepared. The E-glass fibres were impregnated with a ML506 epoxy resin and HA11 hardener. The laminates were fabricated with the hand lay-up technique, and the pre-crack length was produced by positioning a 13 ^m thick Teflon insert at the mid-plane of the plate. In order to produce plates with the desired fibre volume fraction, special pressure tool was applied in order to squeeze out excess resin. Then, the plate was left at room temperature for 24 h to dry. After that, the plate was trimmed with a water jet machine along the longitudinal direction in order to obtain narrow specimens with the desired dimension and initial crack length. After trimming, the nominal length (l) and the nominal width (b) of the DCB specimen were 108 and 25 mm, respectively. The initial crack length (a) was 40 mm. Fibre volume fractions, Vf , measured using the resin burn-off method. Table 1 lists the mechanical property of E-glass/epoxy woven fabrication with Vf = 49% that was used in this research. All of the tests on DCB specimen were Fig. 4. a) DCB specimen, b) typical DCB test, c) crack length measurement during propagation Cohesive Zone Parameters Selection for Mode-I Prediction of Interfacial Delamination 511 Strojniski vestnik - Journal of Mechanical Engineering 61(2015)9, 507-516 completed on a ZWICK electro-mechanical loading frame at room ambient temperature. Fig. 4b shows a typical DCB test. All DCB tests were carried out using displacement control at the crosshead speed of 1 mm/min according to ASTM D5528 standard [19]. Three specimens with a = 40 mm crack length were tested. A load-displacement response was recorded for each specimen during the test. In this work, the crack length was monitored by bonding a strip of paper with the graduations printed on it to the specimen's edge and by taking photos during the experiments with 5 s intervals using a 5 megapixel digital camera. The experimental magnitudes of P-S-a as a function of time were determined. The time of each P-S data point was computed using the applied displacement and the corresponding loading rate. The time for each value of crack length is the one at which the corresponding photo was taken. The photo in Fig. 4c was taken during a test and shows the crack tip, allowing the crack length measurement. These experimental results were then used to verify the adequacy of the three-dimensional finite element scheme utilized to obtain GIc. 3.1.2 Data reduction Method for G Ic In order to calculate the mode-I critical strain energy release rate, there are many commonly used data reduction methods, including compliance calibration method (CCM) that is based on the Irwin-Kies principle, direct beam theory (DBT) and the Corrected Beam Theory. In the present work, the corrected beam theory proposed by de Moura et al. [20] was used. Only one specimen is sufficient to obtain the resistance curve (R-curve) of the specimen, which is the main advantage of the presented method. According to this theory, mode-I critical strain energy release rate can be computed as: GIC = 3PÖ 2b(a + |A| )' (13) where P, a and b are load, crack length and specimen width, respectively. The parameter A is the crack length correction to account the crack tip rotation and deflection. According to the beam theory, the relationship between the compliance and the crack length is expressed by: C1/S = 2(a~ H ) (14) h( Eb)Ui The crack length correction can be obtained using the linear regression of C1/3 versus crack length data. 3.2 Normal Cohesive Strength Measurement The objective of the NCS test is to determine the Mode I cohesive strength (oIc) as an essential parameter for description of the traction separation law of cohesive elements. The fabrication of an NCS specimen is similar to that of a double cantilever beam except for the delaminated area. In the process of NCS specimen fabrication, first a 14-layer woven rectangular plate [0/90] 14 was produced. After the 7th layer of fabrication, a 13 p,m thick Teflon layer was inserted at all sides of the plate so that a 10 mm x 10 mm square area at the middle of plate was released. This area is the cohesive area of NCS specimen, as shown in Fig. 5a. After that, the plate was trimmed to a 80 x 50 rectangular dimension so that the cohesive area was located at the middle of rectangular specimen. The specimen was bonded to fixture surfaces. Prior to bonding, the surfaces of both fixture and specimen, were lightly roughened with the sandpaper on the bonding face and cleaned with acetone. The area of bond is large enough compared to the cohesive area so that debonding would not occur between the specimen and fixture surfaces the testing procedure. Fig. 5a shows the NCS specimen after testing, and Fig. 5b shows a typical NCS test. a; I b) Fig. 5. a) NCS specimen after test, b) NCS test As shown in Fig. 5a, the dimensions of rectangular plates B1 and B2 are 80 mm and 50 mm, respectively, and the width of cohesive the square area (C) is 10 mm. All of the NCS tests were carried out using displacement control at the crosshead speed of 0.5 mm/min until fracture. In this process, three NCS failure tests were completed. In order to prevent slippage during the test, the specimen was accurately balanced and a very little clamping force was required. Decohesion between the fibre and matrix in the NCS cohesive area is the dominant failure mode. Thus, it is obvious that the bulk matrix behaves differently than the thin cohesive layer due to constraint effects induced by the adhesion between the fibre and matrix. 512 Moslemi, M. - Khoshravan, M. Strojniski vestnik - Journal of Mechanical Engineering 61(2015)9, 507-516 As a result, the bulk matrix properties could not be used to determine normal cohesive strength, and this parameter should be determined using an NCS test method. For the purposes of data reduction, all specimens were assumed to have failed instantly. The specimen failure is assumed at its maximum load value. The normal cohesive strength values (oIc) were computed using Eq. (15) as: P — max / (15) where Pmax is the maximum load in which failure occurs and ANCS is the cohesive area of the NCS specimen. 4 RESULTS AND DISCUSSION 4.1 Experimental Results Fig. 6 shows the load-displacement response of three NCS experiments. The measured load is initially negligible, corresponding to the clearance of the fixture. The test was preceded until a maximum load was achieved and followed by a sudden load drop, indicating specimen failure. Displacement [mm] Fig. 6. Load-displacement response of NCS test The mean value of maximum loads of Fig. 6 was considered as Pmax. Using Eq. (15) for data reduction and substituting this value for Pmax, the cohesive strength of composite material computed and is equal to 12.42 MPa. To calculate the experimental resistance curve, the numerical values of P- 8- a parameters were recorded during crack propagation and were used to obtain the critical fracture energies versus crack length. Fig. 7 shows the experimentally obtained R-curve of the material. To simulate the crack propagation using cohesive elements, the mean value of GIc was considered as the fracture toughness of the material. Table 2 lists the cohesive properties of E-glass/epoxy woven composite laminate. Fig. 7. R-curve of E-glass/epoxy specimen Table 1. Mechanical properties of E-glass/epoxy E1 E2 E3 [GPa] [GPa] [GPa] v12 G12 [GPa] Gi3 [GPa] G23 [GPa] 18.43 18.43 3.57 0.15 2.85 2.85 2.85 Table 2. Interfacial property of E-glass/epoxy °Ic TIIc TIIIc [MPa] [MPa] [MPa] GIc [J/m2] GIIc [J/m2] GIIIc [J/m2] n 12.42 22.64 22.64 604 720 720 1.8 4.2 Cohesive Zone Model Results FEM models of each specimen were carried out using the Ply elastic properties of adherent layers that are given in Table 1 and the interfacial properties obtained previously and listed in Table 2. It should be mentioned that when using Eq. (2) for interfacial penalty stiffness, the value of K = 85 x106 N/mm3 is used for all DCB simulations. 4.2.1 Determination of Cohesive Zone Length Cohesive zone length was previously introduced as the distance between the maximum traction and the crack front. Therefore, in order to calculate the distribution of traction in the cohesive layer of model and the corresponding cohesive zone length, a very refined mesh using element size of 0.125 mm was used in the area near the crack tip of the DCB specimen. The distribution of tractions ahead of the crack tip at the delamination onset of the DCB specimen is illustrated in Fig. 8. At the crack initiation point, traction at the crack tip vanished as expected from the cohesive zone theory. According to this analysis, the cohesive zone length of material is 3.76 mm, as shown in Fig. 8. Using the material property that shown in Table 1 and 2 and Eq. (11), the parameter R is computed as 0.27. This value is closest to the Irwin [17] (0.31) model. The cohesive zone length is a material property that has a high order of importance regarding obtaining a successful prediction of delamination onset. This Cohesive Zone Parameters Selection for Mode-I Prediction of Interfacial Delamination 513 Strojniski vestnik - Journal of Mechanical Engineering 61(2015)9, 507-516 parameter was previously introduced as a function of normal cohesive strength by Turon et al. [11]. Thus, using an absolute value for normal cohesive strength, this parameter as a material property is determined. due to the magnitude of tractions ahead of the crack tip. P Distance from crack tip [mm] Fig. 8. Distribution of traction ahead of crack tip 4.2.2 Investigation of Mesh Refinement In order to investigate the effect of mesh refinement in the cohesive zone length on numerical prediction of delamination onset, several DCB specimens were simulated with different sizes of cohesive elements in the cohesive zone length. The predicted load-displacement responses obtained using DCB models are compared to the experimental solution in Fig. 9. In this analysis, cohesive element sizes (in the direction of crack propagation) range from 0.125 mm to 2 mm. The results illustrate that for all mesh sizes a converged solution was obtained but it is necessary to apply a mesh size, le, less than 1 mm to accurately predict delamination initiation. The analysis using coarser meshes significantly overpredicts the strength of the DCB specimen, and the response does not follow the experimental results. Cohesive zone length, lcz, for the material described in Tables 1 and 2 was determined as 3.76 mm. Therefore, for a mesh size smaller than 1 mm, more than three elements would span the cohesive zone length, which is enough for an accurate prediction of the fracture onset. There are several guidelines for the minimum required element in cohesive zone length. For example, Moes and Belytschko [21] proposed using more than 10 elements, while Falk et al. [16] suggested between 2 and 5 elements. In the work of Davila et al. [22], the minimum required element length to predict the delamination in a DCB model was 1 mm, and using more than 3 elements in cohesive zone length of simulated DCB specimen was recommended. The difference in predictions from using a coarse mesh in the modelling of delamination in a DCB specimen is 4 6 Displacement [mm] Fig. 9. Damage simulation using different mesh refinement 4.2.3 Comparison of Damage Models A study also was conducted to investigate the adequacy of the two mentioned traction separation laws used to simulate damage propagation. The objective was to determine how the used models reproduce crack initiation and propagation. Fig. 10 shows the load-displacement of the cohesive zone model and experimental work on a DCB specimen. Fig. 10. Load displacement response of experimental and damage models In the current study, the value of shape parameter a in the PPR model varied from 1 to 4; the result was that increasing the value of the shape parameter a increased the rate at which material loses its stiffness once damage was initiated. In other words, increasing parameter a decreases the fracture process zone effect ahead of the crack tip. In the case of a < 2, there is a gradual fall in the load, but in the case of a > 2 a sudden drop in the load-displacement response is achieved, which means a large number of cohesive elements failed at the same time. For clarity, the results are not shown in this figure. In this study, a 514 Moslemi, M. - Khoshravan, M. Strojniski vestnik - Journal of Mechanical Engineering 61(2015)9, 507-516 value of a = 1.7 was found to be the optimum value for the numerical prediction of damage propagation. As shown in Fig. 10, the modified PPR model was found to be adequate to reproduce the experimentally observed behaviour of the composite specimen, and reproduced approximate smooth crack initiation and propagation while the bilinear model depicted sudden damage propagation. The maximum difference between the experimental and bilinear models is 8.8 % while for PPR it is 2.6 %. This means the modified PPR model accounted fracture process zone which created ahead of crack tip. 5 CONCLUSION A methodology for the delamination characterization of composite laminates under pure Mode I was proposed. • An NCS test has been proposed to compute the Mode-I cohesive strength as a cohesive parameter. • The Mode-I critical strain energy release rate versus the crack length of E-glass/epoxy composite laminate was computed using corrected beam theory for data reduction. • A mixed-mode triangular constitutive relationship between stress (c) and relative displacements (5) of cohesive elements and modified PPR damage model were considered to simulate delamination onset and propagation. • The results of the three-dimensional finite element analysis with cohesive parameters (cIc, GIc) enclosed the adequacy of cohesive parameters. • Accurate damage prediction was achieved using the modified PPR model, and it was considered by the authors to be an accurate model for damage characterization of material. • Modified PPR models accurately described the fracture process zone, which was created ahead of the crack tip as compared to bilinear model. • To ensure the sufficient dissipation of energy, cohesive zone length as a material property was determined. • Numerical analysis with different discretizations of the cohesive zone length showed that numerical predicted responses correlate well with the experimental solutions when at least 3 elements span the cohesive zone length. 6 ACKNOWLEDGMENT This work has been funded by University of Tabriz, and its authors would like to thank the University of Tabriz for the grant. 7 REFERENCES [1] Karpiliski, A. (2006). An Introduction to the diagnosis of the delamination process for glass/epoxy com posites during high-pressure abrasive water-jet cutting. Strojniški vestnik -Journal of Mechanical Engineering, vol. 52, no. 7, p. 532-538. [2] Hasiotis, T., Badogiannis, E., Tsouvalis, N.G. (2011). Application of ultrasonic C-scan techniques for tracing defects in laminated composite materials. Strojniški vestnik - Journal of Mechanical Engineering, vol. 57, no. 3, p. 192-203, D0l:10.5545/sv-jme.2010.170. [3] Brunner, A., Blackman, B., Davies, P. (2008). A status report on delamination resistance testing of polymer-matrix composites. Engineering Fracture Mechanics, vol. 75, no. 9, 2779-2794, D0I:10.1016/j.engfracmech.2007.03.012. [4] Dugdale, D. (1960). Yielding of steel sheets containing slits. Journal of the Mechanics and Physics of Solids, vol. 8, no. 2, p. 100-104, D0I:10.1016/0022-5096(60)90013-2. [5] Barenblatt, G. (1962). The mathematical theory of equilibrium cracks in brittle fracture. Advances in Applied Mechanics, vol. 7, no. 1, p. 55-129, D0I: 10.1016/j.ech.2007.03.012. [6] Hillerborg, A., Modeer, M., Petersson, P.-E. (1976). Analysis of crack formation and crack growth in concrete by means of fracture mechanics and finite elements. Cement and Concrete Research, vol. 6, no. 6, p. 773-781, D0I:10.1016/0008-8846(76)90007-7. [7] Song, S.H., Paulino, G.H., Buttlar, W.G. (2006). Simulation of crack propagation in asphalt concrete using an intrinsic cohesive zone model. Journal of Engineering Mechanics, vol. 132, no. 11, p. 1215-1223, D0I:10.1061/(ASCE)0733-9399(2006)132:11(1215). [8] Zhang, Z.J., Paulino, G.H. (2005). Cohesive zone modelling of dynamic failure in homogeneous and functionally graded materials. International Journal of Plasticity, vol. 21, no. 6, p. 1195-1254, D0I:10.1016/j.ijplas.2004.06.009. [9] Khoshravan, M.R., Moslemi, M. (2014). Investigation on mode III interlaminar fracture of glass/epoxy laminates using a modified split cantilever beam test. Engineering Fracture Mechanics, vol. 127, p. 267-279, D0I:10.1016/j. engfracmech.2014.06.013. [10] Song, S.H., Paulino, G.H., Buttlar, W.G. (2006). A bilinear cohesive zone model tailored for fracture of asphalt concrete considering viscoelastic bulk material. Engineering Fracture Mechanics, vol. 73, no. 18, p. 2829-2848, D0I:10.1016/j. engfracmech.2006.04.030. [11] Turon, A., Davila, C.G., Camanho, P.P., Costa, J. (2007). An engineering solution for mesh size effects in the simulation of delamination using cohesive zone models. Engineering Fracture Mechanics, vol. 74, no. 10, p. 1665-1682, D0I:10.1016/j.engfracmech.2006.08.025. [12] Yuan, H., Li, X. (2014). Effects of the cohesive law on ductile crack propagation simulation by using cohesive zone models. Engineering Fracture Mechanics, vol. 126, p. 1-11, D0I:10.1016/j.engfracmech.2014.04.019. [13] Bhattacharjee, T., Barlingay, M., Tasneem, H., Roan, E., Vemaganti, K. (2013). Cohesive zone modelling of mode I tearing in thin soft materials. Journal of the Mechanical Cohesive Zone Parameters Selection for Mode-I Prediction of Interfacial Delamination 515 Strojniski vestnik - Journal of Mechanical Engineering 61(2015)9, 507-516 Behavior of Biomedical Materials, vol. 28, p. 37-46, DOI:10.1016/j.jmbbm.2013.07.015. [14] Benzeggagh, M., Kenane, M. (1996). Measurement of mixed-mode delamination fracture toughness of unidirectional glass/epoxy composites with mixed-mode bending apparatus. Composites Science and Technology, vol. 56, no. 4, p. 439449, DOI:10.1016/0266-3538(96)00005-X. [15] Park, K., Paulino, G.H., Roesler, J.R. (2009). A unified potential-based cohesive model of mixed-mode fracture. Journal of the Mechanics and Physics of Solids, vol. 57, no. 6, p. 891-908, DOI:10.1016/j.jmps.2008.10.003. [16] Falk, M.L., Needleman, A., Rice, J.R. (2001). A critical evaluation of cohesive zone models of dynamic fracture. Le Journal de Physique IV, 11(PR5), Pr5-43-Pr45-50. [17] Irwin, G. (1997). Plastic zone near a crack and fracture toughness. [18] Abaqus/CAE user manual. Abaqus 6.12 Documentation. [19] ASTM D5528:2007. Standard Test Method for Mode I Interlaminar Fracture Toughness of Unidirectional Fiber-Reinforced Polymer Matrix Composites, ASTM International West Conshohocken. [20] De Moura, M., Campilho, R., Gongalves, J. (2008). Crack equivalent concept applied to the fracture characterization of bonded joints under pure mode I loading. Composites Science and Technology, vol. 68, no. 10-11, p. 2224-2230, DOI:10.1016/j.compscitech.2008.04.003. [21] Moes, N., Belytschko, T. (2002). Extended finite element method for cohesive crack growth. Engineering Fracture Mechanics, vol. 69, no. 7, p. 813-833, DOI:10.1016/S0013-7944(01)00128-X. [22] Davila, C.G., Camanho, P.P., de Moura, M. F. (2001). Mixed-mode decohesion elements for analyses of progressive delamination. Proceedings of the 42nd AIAA/ASME/ASC, DOI:10.2514/6.2001-1486. 516 Moslemi, M. - Khoshravan, M.