Cohesive Zone Parameters Selection for Mode-I Prediction of Interfacial Delamination

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 Cohesive Zone Parameters Selection for Mode-I Prediction of Interfacial Delamination Moslemi, M. – Khoshravan, M. 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


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 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 fibrereinforced 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 (G Ic ) and mode-I cohesive strength (σ Ic ) 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.

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.Fracture simulation of materials using cohesive elements requires substantial experience for determining the mesh requirements and the accurate values of parameters that characterize the tractionseparation law.In this work, in order to simulate Mode-I fracture of E-glass/epoxy woven composite laminate using cohesive elements, two types of traction separation law, including mixed-mode Triangular traction-separation response and modified PPR model, which is potential-based model, were used and compared to experimental results.The objective is to determine an appropriate methodology to predict interlaminar crack growth in composite laminates.As shown in Fig 1, the subscript m referred to mixed mode and subscript c and u referred to critical and ultimate (failure) values, respectively, and i = I, II and III referred to pure modes of fracture.The area under the curves represents the corresponding critical strain energy In order to determine the mixedmode traction-separation law, the properties required to be defined are the three critical strain energy release rates (G ic ), the penalty stiffnesses (K i ) and the cohesive strengths (σ ic ).In the present work, Mode III is assumed to be identical to Mode II, so that the shear strengths in the two orthogonal directions, σ IIc and σ IIIc also critical strain energies G IIc , G IIIc are equal.In the case of the simulation of Mode-I fracture, in order to predict delamination onset, analysis is more sensitive to the parameters σ Ic and G Ic than other interfacial properties.Therefore, in this work, the double cantilever beam (DCB) test and the normal cohesive strength (NCS) test were carried out to determine G Ic and σ Ic , respectively.Other fracture parameters were determined by the calibration of cohesive parameters.In other words, numerical results are fit with experimental results according to the methodology described in [7], [10] and [11].

Damage Initiation Criterion
As shown in Fig. 1, the initial response of the cohesive elements at each damage model is based on linear elastic fracture mechanics (LEFM) and is assumed to be linear until a crack initiation criterion is satisfied.The penalty stiffness, K i , of each traction-separation response law that relates traction to the separation of cohesive elements before crack initiation is defined as below: where i = I, II and III is fracture modes, σ ic and δ ic are the cohesive strength and critical separation of pure modes of fracture, respectively.In the present work, the penalty stiffness for all modes of fracture are considered to be the same, i.e., K I = K II = K III = K.Several methodologies have been proposed to determine the penalty stiffness of cohesive elements.The magnitude of the penalty stiffness must be high enough to avoid interpenetration of the crack surfaces and to avoid artificial compliance from being defined into the model by the cohesive elements.However, a more than enough value can lead to numerical problems.Turon, et al. [11] assumed that whenever the through-the-thickness Young's modulus of the adjacent sub-laminate, E 3 , is small enough compared to K×t, the effective elastic properties of the material will not be affected by the cohesive face.Where t is the thickness of adjacent sub-laminate of composite specimen.Therefore, an equation to calculate the interface stiffness for Mode I is suggested, as mentioned below: where α is a parameter larger than 1.Turon et al. [11] recommended considering α = 50.This value is used for derivation of the penalty stiffness in the analysis presented in this work.After the linear part of the constitutive law, traction at crack tip reaches its maximum value and consequently damage initiates.In order to simulate delamination onset in cohesive elements, there are four well-known crack initiation criteria: maximum stress, quadratic stress, maximum strain and quadratic strain criterion.In this analysis, the quadratic stress criterion that is expressed in Eq. ( 3) was used.Under mixed-mode loading, an interaction between modes must be taken into account.This criterion model readily takes into account the interaction of the traction components in the prediction of damage onset.Moreover, due to the high sensitivity of failure initiation to the strain and displacement, a stress-based criterion gives an accurate failure prediction as compared to other models, such as strain-based or displacement-based criteria.The quadratic stress criterion is defined as: where x, y, z refer to Cartesian coordinate as shown in Fig. 4a.The Macaulay bracket, < >, 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.

Damage Propagation Criterion
After damage initiation, the softening procedure occurs.This procedure is governed by the corresponding traction separation law as below.

Bilinear Traction-Separation Model
As shown in Fig. 1a, in the mixed-mode constitutive law, σ mc and δ 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: where d is a variable that relates damage condition, which has the magnitude d = 0 when the interface 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).

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]: where δ is normal displacement, and A, w, α, δ u are PPR parameters that can be determined by satisfying all boundary conditions, including: where G is the total strain energy release rate.By applying the above-mentioned boundary conditions, it can be expressed as: where δ u can be determined using the following equation as: 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.

Cohesive Zone Length
As shown in Fig. (2), the length of the cohesive zone l cz 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).
where G Ic is the critical strain energy release rate, σ Ic is the normal cohesive strength, and R is a parameter 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.

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: where l e is the length of the cohesive elements in the direction of crack propagation and l cz 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. 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, V f , measured using the resin burn-off method.Table 1 lists the mechanical property of E-glass/epoxy woven fabrication with V f = 49% that was used in this research.All of the tests on DCB specimen were 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-δ-a as a function of time were determined.The time of each P-δ 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 threedimensional finite element scheme utilized to obtain G Ic .

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: where P, a and b are load, crack length and specimen width, respectively.The parameter ∆ 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: The crack length correction can be obtained using the linear regression of C 1/3 versus crack length data.

Normal Cohesive Strength Measurement
The objective of the NCS test is to determine the Mode I cohesive strength (σ Ic ) 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 7 th layer of fabrication, a 13 μm thick Teflon layer was inserted at all sides of the plate so that a 10 mm × 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 × 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.As shown in Fig. 5a, the dimensions of rectangular plates B 1 and B 2 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.
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 (σ Ic ) were computed using Eq. ( 15) as: where P max is the maximum load in which failure occurs and A NCS is the cohesive area of the NCS specimen.

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.The mean value of maximum loads of Fig. 6 was considered as P max .Using Eq. ( 15) for data reduction and substituting this value for P max , 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 - δ - 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 G Ic was considered as the fracture toughness of the material. [GPa] G 13 [GPa] G 23 [GPa] τ IIc [MPa] τ IIIc [MPa] G Ic G IIIc

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 ×10 6 N/mm³ is used for all DCB simulations.

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 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.

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 loaddisplacement 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, l e , 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, l cz , 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, Moës and Belytschko [21] proposed using more than 10 elements, while Falk et al. [16] suggested between 2 and 5 elements.In the work of Dàvila 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 due to the magnitude of tractions ahead of the crack tip.

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.In the current study, the value of shape parameter α in the PPR model varied from 1 to 4; the result was that increasing the value of the shape parameter α increased the rate at which material loses its stiffness once damage was initiated.In other words, increasing parameter α decreases the fracture process zone effect ahead of the crack tip.In the case of α < 2, there is a gradual fall in the load, but in the case of α > 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 value of α = 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.

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 (σ) and relative displacements (δ) 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 (σ Ic , G Ic ) 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.

ACKNOWLEDGMENT
This work has been funded by University of Tabriz, and its authors would like to thank the University of Tabriz for the grant.

Fig. 1 .
Traction-separation law a) triangular model b) modified PPR model

Fig. 3 .
Fig. 3. Deformed shape of simulated DCB with 3D elements 3 EXPERIMENTAL APPROACHES 3.1 Mode I Fracture Toughness Measurement 3.1.1The DCB Test Fig. 4. a) DCB specimen, b) typical DCB test, c) crack length measurement during propagation 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].

Fig. 8 .
Fig. 8. Distribution of traction ahead of crack tip

Fig. 10 .
Fig. 10.Load displacement response of experimental and damage models

Table 2
lists the cohesive properties of E-glass/epoxy woven composite laminate.

Table 1 .
Mechanical properties of E-glass/epoxy

Table 2 .
Interfacial property of E-glass/epoxy