Finite Element Simulation of Interlaminar and Intralaminar Damage in Laminated Composite Plates Subjected to Impact

Here In this study, the composite laminates subjected to transverse impact with consideration interlaminar and intralaminar damage based on Cohesive Zone Model (CZM) and Progressive Damage Model (PDM) are investigated by numerical analysis using ABAQUS commercial finite element code. The delamination in stacking ply with the same fiber orientation is considered as interlaminar damage and the delamination in an inner layer of any cluster is ignored. Hashin criterion is used for intralaminar damage initiation and evolution without using any subroutine. First, the appropriate procedure for delamination on composite specimen was suggested based on CZM approach in double cantilever beam to verify the intralaminar damage simulation. Then by considering several case studies with different impact energies, the results of present simulation is verified with the relevant and available experimental results and numerical references in the existing literature. According to the available experimental results the present simulation results are more acceptable and accurate than the results of similar numerical works, especially in higher impactor velocity.


INTRODUCTION
Composite structures are widely used in many industrial applications like aerospace and defence industry due to their inherently high specific mechanical properties. In many situations, these composite structures are subjected by high and low velocity impact. Composite structures are very sensitive to non-visual damage that strongly influence their residual load bearing capability. Lack of knowledge of the impact effects on composite structures is a factor in limiting the use of composite materials (Abrate, 1998). To understand the responses of composite structures under high velocity projectile impacts, various experimental and numerical studies have been conducted. Cantwell and Morton (1989) examined the initiation and development of damage in composite structures with a series of low and high velocity impact tests. Guoqi et al., (1992) investigated experimentally the response of woven Kevlar/polyester laminates of varying thicknesses to quasi-static and dynamic penetration by projectiles. Cheng et al., (2003) developed a model based on hydrodynamic finite element code for high velocity impact on thick composites. This model was based on a continuum approach which was built on the basis of an orthotropic constitutive behaviour with stress-based failure criteria and a simplified degradation model of the failure of composites. Silva et al., (2005) investigated experimental and numerical simulation of ballistic impact on Composite structures made of Kevlar-29. Cerioni (2009) presented a numerical simulation of delamination in composite materials under static and fatigue loading by cohesive zone models. Zhao et al., (2010) experimentally reported the failure modes of T300/epoxy composite laminates at different impactor velocities of 10-300 m/s on the different stacking sequence. Khalili et al., (2011) proposed a verified finite element model using ABAQUS finite element for composite laminates and shell structures subjected to low-velocity impact. Gonzalez (2011) investigated the damage induced in composite plates under drop-weight impact loading by analytical description and experimental test plan for assessment of the virtual test performance by finite element simulation. Ramadhan et al., (2013) investigated the high velocity impact response of composite laminated plates, both experimentally and numerically.
Simulation of damage induced in composite laminates has been performed with macro, meso and micro-scale of modelling. Investigating the overall response of the composite laminate with treating the composite as fully homogenized is called marco-scale analysis. In the meso-scale approach the composite treated as effective anisotropic

SIMULATION OF DELAMINATION ON DCB COMPOSITE SPECIMEN BASED ON CZM
In this section, the valid experimental and numerical reference are used to verify the proposed simulation approach. The results of the present simulation are compared with the experimental and numerical results reported by Cerioni (2009). Consequently, the appropriate procedure was suggested for delamination on composite specimen based on CZM approach. An unidirectional 24-ply Graphite/Epoxy composite laminates are considered with the dimensions of 140×20 mm, total thickness of 3.24 mm, laminate configuration [011/-5/5/011] and initial crack length. The DCB (Double Cantilever Beam) test is one of the most common tests used to evaluate the mode I interlaminar fracture toughness in a composite laminate. In fact, the test considers a composite beam with a cohesive layer of thickness t = 0.02 mm and an initial delamination crack of 50 mm length. As shown in the Figure 1 the initial delamination is forced to open by applying a displacement that pull the two beams of the specimen away from each other. The material properties given for composite plies of DCB specimen and cohesive interface are indicated in Table 1 and 2 (Cerioni, 2009).   As indicated in Figure 2, the cohesive zone model, as shown combines an initially linear elastic behaviour with strength-based failure criterion to predict the damage initiation and a fracture mechanics-based criterion to determine the damage evolution. The elastic behaviour is written in terms of an elastic constitutive matrix that connect cohesive surface tractions (t) to displacement (δ) at the interface. The initial stiffness of CZM (Figure 2, E0) defined as traction divided by separation, (units N/m 3 ) should be chosen as high as possible to provide a reasonable stiffness. But from a numerical perspective, it cannot be infinitely large; otherwise, it leads to numerical ill-conditioning (Turon et al., 2007;Dssault System's Simulia Corp, 2016) Many researchers have suggested various guidelines for selecting the stiffness of the interface. For instance, Camanho and Matthews (1999) obtained predictions for graphite/epoxy specimens by using a value of 10 6 N/mm 3 . Zou et al., (2002) proposed a value for the interface stiffness between 10 4 and 10 7 times of the interfacial strength value per Latin American Journal of Solids and Structures, 2018, 15(6), e90 4/20 unit length. Turon et al., (2007) proposed equation (1) in order to estimate the stiffness of the interface and declared that the equation is more appropriate than the other works. In this equation, E3 is the elastic modulus through-the-thickness of the composite materials t is the laminate thickness adjoining to the cohesive interface, and α is an increasing parameter normally taken about 50.
As Shown in Figure 2, the point t0 refers to the beginning of materials response degradation. The t0 represents the peak values of the traction that correspond to the δ0 when the deformation is purely normal to the interface. The process of degradation begins when the traction or separation satisfies certain damage initiation criteria. Several damage initiation criteria such as maximum nominal stress (MAXE Damage) and quadratic nominal stress (QUADS Damage) are available in the ABAQUS (Dssault System's Simulia Corp, 2016). Each damage initiation criterion also has an output variable associated with it in order to indicate whether the criterion is met. Perillo et al., (2012) have reached a better agreement between experimental results and the finite element results by using the QUADS initiation criterion. The Quadratic nominal stress criterion is assumed to initiate when a quadratic interaction function reaches a value of one. This criterion can be represented as equation (2). In three-dimensional problems, the index n refers to normal traction and the index s and t refer to two shear traction. The Macaulay bracket (‹ ›) are used to signify that a pure compressive deformation or stress state does not initiate damage. Under single-mode loading, interface damage initiates when the traction reaches the maximum nominal interfacial strength (t 0 ) (Dssault System's Simulia Corp, 2016). According to the strength-based failure criterion, the damage evolution would be initiated once the damage initiation condition has been met. Several constitutive models that have been proposed by the researchers for the damage evolution schemes, linear softening behaviour is usually implemented (Camanho et al., 2003). The damage evolution law describes the rate at which the material stiffness is degraded once the corresponding initiation criterion is reached. A scalar damage variable, D, represents the overall damage in the material. It initially has a value of 0 and then monotonically evolves from 0 to 1 upon further loading. The stress components of the traction-separation model are affected by the damage. According to equation (3), t is the stress components predicted by the elastic traction-separation behaviour for the current strains without damage (Dssault System's Simulia Corp, 2016). As shown in Figure 2, a value of 1 of the SDEG (Overall value of the scalar damage variable D), indicates that the initiation criterion has been met at the end of the damage evolution area which corresponds to the δm. Camanho et al., (2003) proposed equation (4) in order to define the effective displacement (δm) to describe the evolution of damage under a combination of normal and shear deformation across the interface.
Damage evolution can be defined based on the energy that is dissipated as a result of the damage process. The area under the traction-separation relation (the area that shaded in Figure 2) is equal to the fracture toughness energy, Gc. Although the scalar damage variable (D) is derived from this parameter, the fracture energy (Gc) is the most important parameter which can be determined by means of some standard experimental tests.

Solution method
The simulation of delamination based on CZM approach, is possible to be solved by an explicit or implicit algorithm available in ABAQUS. The ABAQUS/Standard general-purpose solver uses a traditional implicit integration scheme to solve finite element analyses. On the other hand, the ABAQUS/Explicit which uses an explicit centraldifference time integration rule applies an explicit integration scheme to solve highly nonlinear transient dynamic and quasi-static analyses. The simulation of delamination based on CZM approach is involves softening in the material response, and lead to convergence difficulties in an implicit solution procedure in ABAQUS/Standard. Convergence difficulties may also occur during unstable crack propagation when the available energy is higher than the fracture toughness of the material (Dssault System's Simulia Corp, 2016). Khoramishad et al. (2010) uses standard solver of ABAQUS to simulate the single lap joints according to the CZM approach. Cerioni (2009) offers explicit solver in order to obtain the finite element solution for avoiding convergence problems and numerical instabilities. In the earlier work (Khalili et al., 2011), it was found that in an explicit scheme, the analysis cost rises only linearly with problem size, whereas the cost of solving the nonlinear equations associated with implicit integration rises more rapidly than linearly with problem size. Therefore, ABAQUS/Explicit is attractive for very large problems.

Type of elements
Based on plane strain assumption, the simulation of delamination on DCB composite could be analysis with twodimensional elements. Mi et al., (1998) applied a three-dimensional analysed and indicated that the 2D plane strain formulation gave a very reasonable approximation versus three-dimensional analysis. While Alfano and Crisfield (2001) strictly declared that a 3D analysis would be needed since the delamination front was not precisely a straight line.
Conventional shell, continuum shell and solid elements are available options in ABAQUS for modelling the composite structures. While cohesive behaviour is used in conjunction with stacked conventional shell elements, the specialized contact formulation may lead to approximate normal contact forces. This may induce approximate transverse shear behaviour in the stacked shells which affect the bending behaviour of the stack. Continuum shells should be used instead of conventional shells in such modelling scenarios (Dssault System's Simulia Corp, 2016).

Cohesive surfaces versus cohesive elements
Similar to cohesive surfaces, cohesive elements could be used in CZM approach simulation in ABAQUS. The formulation used for surface-based cohesive, is very similar to that of cohesive elements with traction-separation response. For cohesive surfaces, the cohesive constraint is enforced at each slave node while in cohesive elements, the cohesive constraints are calculated at the material points. Hence for cohesive surfaces, refining the slave surface as compared to the master surface will likely lead to improved constraint satisfaction and more accurate results (Dssault System's Simulia Corp, 2016). Turon et al., (2007) used cohesive elements in CZM approach, but others like Gonzalez (2011), used both of cohesive elements and cohesive surfaces in CZM approach and concluded that the use of cohesive surfaces resulted in improving the runtime of the analysis. On the other hand, it can yield to large runtime analysis because contact algorithms are computationally heavy especially when the connected mesh are highly refined.

Verification and discussion on simulation method
This section briefly describes the finite element (FE) simulation of the delamination on DCB composite specimen based on CZM approach. In this study, the FE simulation, has been done with the finite element package ABAQUS version 2016 and is performed on a computer with 3.2 GHz Intel ® Corei7, with 8.00 GB RAM.
According to Table 1, the material modelling of composite ply in the present FE simulation has been done by linear elastic with engineering constant option that is available in ABAQUS. The CZM approach is used to simulate the interface failure modelling of cohesive layer between two beams of the specimen. According to Table 2, the magnitude of the initial stiffness of the cohesive element (Figure 2, E0) has been imported to the ABAQUS and it does not need to be divided by the cohesive thickness. For defining the quadratic nominal stress (QUADS Damage) criterion, the S and T which are taken from Table 2 are considered as the maximum nominal interfacial strength (t 0 ) in the two directions (normal and shear). Also according to Table 2, the fracture energies GIC, and GIIC are imported to ABAQUS to define the modes I and II for simulation of the interlaminar fracture toughness in cohesive interface.
The 4 nodes, two-dimensional continuum plane strain elements (CPE4R) and the 4 nodes, two-dimensional cohesive elements (COH2D4) are used in the FE simulation of the composite specimens and the cohesive interface. A very refined mesh using an element length of 0.25 mm was used for both solid and cohesive elements and is constant for all the models (6720 elements of type CPE4R and 360 elements of type COH2D4). This element length is recommended by Cerioni (2009) and as shown in Figure 3, only one cohesive element in the cohesive interface was put. Here, a double cantilever beam has been simulated under displacement control. As shown in Figure 4, the movement of the opening edge in the vertical direction of the composite specimens is applied by the displacement control and is restrained in horizontal direction. Also, the last node of the bottom of the lower composite specimen has been restrained in vertical direction. In this case study, the present simulation results gained by standard solver of ABAQUS are compared with the available experimental and numerical results by the Cerioni (2009). As shown in Figure 5, there is a good correspondence between the elastic branches of the three curves. Also, the prediction of the onset of the delamination looks adequate. But some mismatch was observed between the present simulation and the results of reference (Cerioni, 2009) in onset of the delamination. Furthermore, there is a difference in the propagation branch of the simulation model stiffness which has been decreased more quickly in respect of the experimental results. In addition, better correspondence was observed with experimental results in the present simulation versus simulation of Cerioni (2009). In Figure 6, the present simulation results gained by the explicit solver of ABAQUS are compared with the available experimental and numerical results by the Cerioni (2009). As indicated in Figure 6, the mismatch in onset of delamination gained by standard solver of ABAQUS ( Figure 5) is eliminated. Also better results correspondence was observed in explicit solver than standard solver. Figure 6: Comparison of contact force variation as a function of displacement for delamination on DCB composite specimen between the present FE simulation with explicit solver and the results of Cerioni (2009) As mentioned earlier in section 2.3, a 3D analysis would be needed since the delamination front is not precisely a straight line (Alfano and Crisfield, 2001). In this study the composite specimen uses the SC8R element for investigation the effect of 3D modelling in delamination. According to the earlier work (Khalili et al., 2011), SC8R elements should be used in thick plates and shells. The SC8R or S4R elements are an eight-node continuum shell and four-node conventional shell elements with reduced integration, respectively.
In case of modelling, cohesive behaviour is possible to be used in conjunction with stacked conventional shell elements. Depending on the load case, the specialized contact formulation may lead to approximate normal contact forces, which in turn may induce the approximate transverse shear behaviour in the stacked shells and consequently affect the bending behaviour of the stack. For increasing accuracy in such modelling scenarios, continuum shells should be used instead of conventional shells (Dssault System's Simulia Corp, 2016).
As mentioned earlier in section 2.4, cohesive surfaces like cohesive elements could be used in CZM approach simulation. Also unlike to 2D modelling, the 3D modelling leads to large runtime analysis. Therefore the cohesive surfaces should be uses in 3D modelling in order to improve the analysis runtime.
In the cohesive surface method, the initial stiffness magnitude of the cohesive surface that should be imported to the ABAQUS is different from the magnitude for the cohesive element. In addition it should be divided by the cohesive thickness and does not need to multiply to 50, according to the equation (1).
As indicated in Figure 7, the force variation as a function of displacement for delamination on DCB composite specimen by using SC8R elements (3D modelling) is similar to Figure 6, in which explicit solver is also used. But the fluctuation in the force history of the 3D modelling versus 2D modeling is considerable. The non-symmetric pattern of the displacement contour which is shown in Figure 8, leads to the fluctuation of force history. Therefore unlike to the 2D modeling, the fluctuation of force history is considerable in the 3D modelling because the delamination is not precisely a straight line.

SIMULATION OF IMPACT ON LAMINATED COMPOSITE PLATES BASED ON CZM AND PDM
By applying the appropriate procedure for delamination based on CZM method mentioned in section 2; here, the impact FE simulation process of composite plates is briefly described regarding CZM and PDM approach. To verify the proposed simulation approach, the results of the present simulation are compared with the valid experimental and numerical results reported by Gonzalez (2011). And finally, the appropriate procedure for impact on composite laminated plates is suggested by regarding intralaminar and interlaminar damage. A plate is considered with the dimensions of 150×100 mm 2 made of Hexply AS4/8552 carbon-epoxy unidirectional pre-preg composite with laminate configuration [454/04/-454/904]s. In Figure 9 (Gonzalez, 2011), the test specimen and fixture base with details of support area and clamping points of the specimen are indicated. The material failure properties assumed for AS4/8552 plies and cohesive interface are demonstrated in Tables 3 and 4 (Gonzalez, 2011). It should be mentioned that the ʋ12 is 0.35, the mass density of the plate is 1.590 kg/m 3 , and the ply thickness is 0.18125 mm. The plate is subjected to transverse impact by a steel sphere impactor of 16 mm diameter with the mass of 5 kg.

Material failure modelling
In this study, the material failure modelling is based on the progressive damage model that is used mainly to obtain the overall response of the composite laminates. A typical PDM consists of three steps: linear elastic stress analysis, failure analysis, and property-degrading material state variable. There are many material models that can be used to predict the behaviour of composite laminates. Progressive damage model in ABAQUS is assumed to be linearly elastic. It is intended that the model predicts the behaviour of composite laminates for damage which can be initiated without a large amount of plastic deformation. The initiation criteria are used to predict the onset of damage and the damage evolution law is based on the energy dissipated during the damage process and linear material softening (Dssault System's Simulia Corp, 2016).
The ABAQUS anisotropic damage model is based on the work of Hashin and Rotem (1973), Hashin (1981) and Camanho et al., (2003). In ABAQUS, the damage initiation criteria for fiber composite laminates are based on Hashin's theory. There are four different damage initiation criteria in ABAQUS (Dssault System's Simulia Corp, 2016) as follows: • fiber rupture in tension • fiber buckling and kinking in compression • matrix cracking under transverse tension and shearing • matrix crushing under transverse compression and shearing The equations for these damage initiation criteria are listed in the following equations (5-8). Tensile or compressive damage can be initiated in the fiber or the matrix when the respective F equals one.
In the above equations, S and T S denotes the longitudinal tensile strength, the longitudinal compressive strength, the transverse tensile strength, the transverse compressive strength, the longitudinal shear and transverse shear strength respectfully. Also,  is a coefficient that determines the contribution of the shear stress to the fiber tensile initiation criterion. In addition 11 12

22
, and    are components of the effective stress tensor which is used to evaluate the initiation criteria. The initiation criteria presented above can be specialized to obtain the model proposed in Hashin and Rotem (1973)  Prior to damage initiation, the material in ABAQUS Progressive damage modelling is linearly elastic with the stiffness matrix of a plane stress orthotropic material. Thereafter, the effect of damage is taken into account by reducing the value of stiffness coefficients. If any of the damage variables become one, stress can no longer be supported in the respective direction because the stiffness goes to zero. The response of the material is computed by the equation (9), where  is the strain and d C is the damaged elasticity matrix (Dssault System's Simulia Corp, 2016).
In a damage evolution model in ABAQUS PDM, it's required to predict the material response in each of the four failure modes as loading increases or unloading. As it's indicated in Figure 10, these should be beyond the point of initial damage as linear damage evolution. Prior to damage initiation the positive slope of the stress-displacement curve corresponds to linear elastic material behavior. Whereas the negative slope is achieved by evolution of the respective damage variables after damage initiation. To relieve mesh dependency during material softening, ABAQUS introduces a characteristic length into the formulation. So the damage evolution is expressed as a stress-displacement relation. In For each failure mode, the energy dissipated should be specified due to failure ( c G ), which corresponds to the area of the triangle OAC in Figure 10. The values of t eq  for the various modes depend on the respective c G values.
Thereafter the damage initiated, unloading from a partially damaged state such as point B in Figure 10, occurs along a linear path toward the origin with degraded stiffness and strength due to damage induced.

Solution method
As mentioned in the earlier work (Khalili et al., 2011), the simulation of impact on composite laminates, is possible to be solved by both explicit and implicit solver of ABAQUS, but the ABAQUS/Explicit is the appropriate choice to solve highly nonlinear transient dynamic and progressive damage modelling (Dssault System's Simulia Corp, 2016). Also as it is mentioned in section 2.2, the ABAQUS/Explicit is attractive for very large problems.

Contact modelling
Hard contact law is one of the contact laws that is available in ABAQUS. In this contact method, the contact constraint is applied when the clearance between two surfaces becomes zero. There is no limit in the contact formulation on the magnitude of contact pressure which can be transmitted between the surfaces. The surfaces are separated when the contact pressure between them becomes zero or negative, then the constraint is removed. The contact force is a function of the penetration distance. And it is applied to the slave nodes to oppose the penetration while equal and opposite forces act on the master surface at the penetration point. When using hard contact, it is still possible to make master surface penetrate into the slave surface by pure master-slave algorithm as indicated in Figure 11 (Dssault System's Simulia Corp, 2016). Both the pure master-slave and the balanced master-slave contact algorithms are available in ABAQUS/Explicit. By default, ABAQUS/Explicit will decide which algorithm to be used for any given contact pair based on the nature of the two surfaces forming the contact pair and also applying kinematic or penalty enforcement of contact constraints. ABAQUS/Explicit uses the pure master-slave in kinematic contact algorithm as a rigid surface contacts a deformable surface. If the penalty contact algorithm is specified, ABAQUS/Explicit chooses balanced master-slave weighting as two element-based surfaces contact each other. Balanced master-slave weighting means that the corrections produced by both sets of contact calculations are weighted equally. The penalty contact algorithm for balanced master-slave contact surfaces computes contact forces that are linear combinations of the calculated pure master-slave forces. One set of forces is calculated by considering one surface as the master surface, and the other forces are calculated by considering the same surface as the slave surface (Dssault System's Simulia Corp, 2016). As mentioned in section 2.4, it is possible to use both cohesive element and cohesive surfaces for defining cohesive behaviour between tow laminates of composite plate. But in the 3D modelling, cohesive surfaces are preferable due to their low-cost analysis.

Type of elements and mesh convergence study
Conventional shell, continuum shell and solid elements are available options in the modelling of the composite laminates in ABAQUS. The choice of an element type used for the plates or shells as target structures depends on their side to thickness ratio as well as the impactor velocity. Increasing the thickness of target and impactor velocity lead to increase in the shear deformation.
In order to obtain the optimum number of elements, an analysis of the number of elements sensitivity is performed. In Figure 12 mesh pattern of the circular area surrounding the contact region and the convergence study of the number of elements are shown. This study indicated that the element number of 15 in edges of circular impact area is the optimum number of elements.

Verification and discussion on simulation method
Similar to the section 2 of this study, the ABAQUS version 2016 is used for simulation of the impact on laminated composite plates. The simulation is performed on the same computer mentioned in the previous section.
The material modelling of composite ply has been done by linear elastic with lamina option that is available in ABAQUS. Due to the plane stress condition that is used for SC8R elements, only the values of 11 E , 22 E , ʋ12, 12 G , 13 G and 23 G are required to define an orthotropic material (Dssault System's Simulia Corp, 2016). Therefore based on G is equal to 12 G and 23 G is equal to 12 2 3 G .
As mentioned earlier in section 3.1, Hashin's damage initiation criteria parameters are required in the PDM approach that is used for intralaminar failure modelling of laminated composite. According to the Table 3, Table 3, as fracture energy in the fiber and matrix in tension and compression mode. According to Table 4, the interlaminar failure properties of interfaces are defined. Also the procedure that is used in this section is similar to section 2.5 for cohesive zone modelling. In order to modelling the interlaminar failure of interfaces in the present laminate configuration [454/04/-454/904]s the clustering ply was chosen. Therefore the delamination in stacking ply with the same fiber orientation (clustering) is considered and the interlaminar failure in inner layer of any cluster is ignored. As indicated in Figure 13, seven cohesive surfaces duo to low-cost analysis are modelled according to the CZM procedure by using the properties of Table 4. Based on the preferences of explicit solver mentioned earlier in section 3.2, ABAQUS/Explicit is used in this study for simulating the impact phenomena. The hard contact law is chosen for being applied as the clearance between two surfaces becomes zero. The impactor and the target are set as the master and slave surfaces, respectively with penalty contact algorithm of balanced master-slave contact surfaces. The impactor was modelled as deformable body with element of C3D8R and the laminated composite are meshed using SC8R elements. As mentioned earlier in section 3.4, the optimum number of elements (12224 elements of type SC8R and 4025 elements of type C3D8R) are obtained based on the mesh convergence study, According to section 3, three case studies with different level of impact energies were chosen to verify the present FE simulation procedure. Therefore, impactor with initial velocities of 3.93, 3.38 and 2.78 m/s were simulated.
The results of impact with impactor velocity of 2.78 m/s and kinetic energy of 19.3 J were compared with experimental and numerical results those presented by Gonzalez (2011). The comparison of the results for the contact force history was indicated in Figure 14. As it is observed, there is appropriate correspondence between the results. Also unlike to the simulation results that presented by Gonzalez the fluctuations of contact force are decreased in the present FE simulation. In Figure 15, the FE simulation and experimental variation of energy dissipation of impactor is demonstrated as a function of time for laminated composite plate. The results indicated that the energy absorption by laminated composite plate of the present simulation has less discrepancies than simulation results presented by Gonzalez. In addition, in Figure 16, the comparison of contact force is indicated as a function of impactor displacement for a composite plate under impact kinetic energy of 19.3 J and appropriate correspondence is observed. The present results for the contact force history with impactor velocity of 3.38 m/s and kinetic energy of 28.6 J is shown in Figure 17, and compared with experimental and numerical results of Gonzalez (2011). It should be mentioned Gonzalez (2011) has only reported the contact force history for this level of impact energy. The comparison of the results for the contact force history, energy absorption and contact force as a function of the impactor displacement is shown in Figure 18, Figure 19 and Figure 20 respectively in case of impactor velocity is 3.93 m/s and kinetic energy is 38.6 J. Here, appropriated correspondence is observed. The comparison of contact force history of three case studies (Figure 14, Figure 17 and Figure 18) is indicated that unlike to the present simulation, increasing impactor velocity leads to increasing the discrepancy between simulation and experimental results of Gonzalez (2011). Therefore the validity of the present impact simulation procedure is not dependent on the impactor velocity. In addition as indicated in Figure 18, the maximum contact force that reported by Gonzalez experimentally (point A) is equal to 8.66 kN. The discrepancy in contact force (point A) is equal to 13.43% and 25.24% for the present and Gonzalez simulation respectively. Therefore, the accuracy of the present impact simulation is prefer than the Gonzalez simulation especially in higher impactor velocity.

CONCLUSION
In this study, a reliable, low-cost and simple method is provided based on ABAQUS/Explicit package by considering damage in laminated composite plates under transverse impact loadings. The best procedure is proposed which can serve as benchmark method in damage modeling of composite structures under high velocity impact for future investigations. First, the simulation of delamination on DCB composite specimen is verified based on CZM approach. In addition, materials model, solution method, element type and method of cohesive definition are considered. It is observed that, less CPU run-time is required for simulation of delamination problems 2D modelling. Furthermore, cohesive surface rather than cohesive element should be used in 3D modelling for improving the analysis runtime. Contrary to cohesive surface, there is no need to divide the initial stiffness of the cohesive element -which should be imported to the ABAQUS-by the cohesive thickness It can be declared that, explicit solver of ABAQUS is the appropriate choice for modelling progressive damage and cohesive zone. In addition ABAQUS/Explicit is capable to solve highly nonlinear transient dynamic and it is attractive for very large problems.
In this study, the simulation of impact on laminated composite plates based on CZM and PDM are verified by the experimental and numerical results available in the literature.
By considering damage evolution behaviours of matrix and fiber cracking and interface delamination in three case studies with different levels of impact energies, our simulation results have an appropriate correspondence with the results of similar works especially in the aspect of force-time, force-displacement and energy time histories curves. According to the simulation results, the delamination in clustering ply is significant but the interlaminar failure in an inner layer of any cluster could be ignored.
Like to Gonzalez (2011), we have used ABAQUS finite element simulation; unlike to him, our simulation results are more accurate -about 12 percent better correspondence in maximum contact force-than his simulation results especially in higher impactor velocity. On the other hand, our results correspond to his experimental results appropriately. Finally, the proposed method can serve as a benchmark for simple impact simulation of composite structures based on CZM and PDM in the future investigations, such as optimization study and engineering application of composite laminates under impact.