Numerical comparison between Hashin and Chang-Chang failure criteria in terms of inter-laminar damage behavior of laminated composite

In this paper, a general three-dimensional finite element model composed of eight-node three-dimension cohesive elements and eight-node solid elements with reduced integration for low-velocity impact analysis of composite materials is established. Firstly, the continuum damage mechanics is applied to simulate the initiation and evolution of the intra-laminar damage. Based on cohesive zone model, the cohesive elements are inserted between layers to predict the inter-laminar damage. Three failure criteria, including 2D Hashin, 3D Hashin and Chang-Chang criteria, are coded in VUMAT and are implemented in ABAQUS/Explicit. The numerical results of energy time curves, force time curves, force displacement curves and damage distribution under three impact energies (7.35 J, 11.03 J and 14.7 J) are in good agreement with previously published data in literature, which indicates that the finite element model is suitable for studying the mechanical response and damage distribution of composites laminates subjected to low-velocity impact. Moreover, the influence of stacking sequence and friction coefficient on mechanical response and damage distribution is analyzed. It is concluded that the composite laminate with a stacking sequence of [45°/0°/−45°/90°]S can reduce the area of damage region compared to [0°/90°]2S because the ±45° layer can improve the shear resistance of composite laminate. Also, the computation accuracy will be the best when the friction coefficient is adopted between 0.5 and 0.7.


Introduction
Due to the excellent characteristics of high specific strength, high specific modulus and low specific gravity [1], composite materials have been widely implemented in aerospace fields. However, the poor impact resistance of composite materials makes the barely visible impact damage (BVID) [2,3] become a critical issue in the design and application of composites. The finite element method (FEM) has become a popular research method, which has obvious advantages in the research of mechanical behavior and damage distribution [4,5]. Therefore, it is necessary to establish a finite element model to predict the low-velocity impact mechanical behavior of composite materials and damage distribution [6], which can provide reference suggestions for the impactresistant design of composite materials.
In numerical analysis, the intra-laminar damage can be simulated by the continuum damage mechanics (CDM), which consist of damage initiation criterion and damage evolution law of composite materials. The damage initiation and damage modes of composite materials are judged by the failure criteria. Among these well developed failure criteria, Hashin [7,8] and Chang-Chang [9,10] criteria are the most popular. After years of development and improvement, the Hashin criterion has expanded from 2D to 3D. Compared with 2D Hashin criterion, the 3D Hashin criterion considers the influence of the through-thickness normal stress component s 33 on matrix damage. The Chang-Chang criterion is obviously different from the Hashin criterion in the expression of the matrix damage criterion. Once the failure criterion is satisfied, the corresponding damage

Methodology
General speaking, the methodology of studying the low-velocity impact of composite materials includes two aspects: (1) Appropriate damage constitutive model to simulate the damage modes and mechanical behavior of composite materials; (2) Reasonable finite element model to ensure calculation accuracy.

Damage models of composite materials
The numerical analysis of composite under low velocity impact includes two aspects: (1) Intra-laminar damage, including appropriate failure criterion to predicted damage modes and effective damage evolution model to describe the composite stiffness degradation; (2) Inter-laminar damage, including suitable method to predict the initiation and evolution of delamination [22].

Failure criteria of composite materials
The Hashin criteria is one of the most popular failure criteria which is widely implemented in predicting the initiation of composite materials and the application field of Hashin criterion has expanded from 2D to 3D. The main difference between 2D Hashin and 3D Hsahin is that the through-thickness normal stress component s 33 is taken into consideration in 3D Hashin, which is more in line with actual stress distribution [23][24][25]. The Chang-Chang criterion is also implemented in many studies to predict the low-velocity impact of composite materials, and its matrix damage criterion is obviously different from the Hashin criterion. However, few studies have systematically compared the effects of three failure criteria on the mechanical response and damage distribution of composite laminates under low-velocity impact. The expressions of the three failure criteria are as follows: (1) 2D Hashin failure criterion [7]: Fiber tensile failure s  0 :  Where s s s , , 11 22 33 are the effective stress tensor in corresponding direction; s s s , , 12 13 23 are shear stress; X X Y Y , , .
T C T C are fiber tensile strength, fiber compressive strength, matrix tensile strength and matrix compressive strength, respectively; S S S , , 12 13 23 are shear strengths.

Damage evolution model
Once a failure criterion is satisfied, the damage evolution model is needed to be defined to describe the progressive damage of composites. The constitutive relationship of the degradation model based on energy evolution is shown in figure 1 and the expression is shown in equation (13). Where the e 11 is the strain calculated by each increment step; X T and X C are the fiber tensile strength and fiber compressive strength of the unidirectional composite laminate in fiber direction. G ft and G fc are fiber tensile fracture energy and fiber compressive fracture energy. E 11 is the elastic modulus in fiber direction. l c is the characteristic length of elements, which is equal to the cube root of the element volume.

Matrix failure modes
It is found that the matrix failure criterion of 2D Hashin criterion and Chang-Chang criterion is only related to the s . 22 Therefore, the e 22 is taken into the calculation of matrix damage variables of 2D Hashin criterion and Chang-Chang criterion. However, the 3D Hashin's matrix failure criterion is judged by s s + .
Where G mt and G mc are matrix tensile fracture energy and matrix compressive fracture energy. E 22 is the elastic modulus in transverse direction. The symbol <> represent the Macaulay operator.

Continuum damage model
The damage variable of each damage modes is needed to be substituted into the composite material constitutive for further calculation to simulate damage initiation and stiffness degradation of composite material. The degraded compliance matrix S d is expressed as: Where the d f is the fiber damage variable and d m is the matrix damage variable. They can be calculated by equation (24). The corresponding degraded stiffness matrix C d is shown as:  21 23 32 13 31 21 32 13 ft fc mt mc represent the fiber tensile damage variable, fiber compressive damage variable, matrix tensile damage variable and matrix compressive damage variable, respectively. S S , mt mc are introduced coefficients to control the shear stiffness loss. In present study, S mt and S mc are set as 0.9 and 0.5, respectively.

Inter-laminar damage
The inter-laminar damage is also called delamination. In this paper, the cohesive elements based on CZM are inserted between plies to simulated the distribution of delamination. The damage initiation is judged by quadratic failure criterion, as shown in equation (25). The damage evolution is controlled by Benzeggagh and Kenane (B-K) criterion [28,29], as shown in equation (26).
C are the total, normal and shear critical fracture energy, respectively. G S is the dissipated energy in the out-of-plane direction;G T is the total dissipated energy in all three directions; h is the relevant material coefficient in the B-K criterion and is taken 1.45 in this paper.

Numerical model of low velocity dynamic impact
The low-velocity impact is a strong nonlinear contact process that occurs in a short period of time. Therefore, establishing an appropriate numerical model is fairly necessary, which can not only ensure the convergence of the numerical analysis, but also reduce the computation time.

Finite element model and boundary conditions
According to the geometric dimensions provided by Shi [21], the low-velocity impact three-dimensional numerical analysis model of composite materials is established, as shown in figure 2. The diameter of the composite circular plate is 75 mm. The thickness of the composite laminate with a stacking sequence of [0/90] 2S is 2 mm. The x, y, z degrees of freedom of the edge are constrained to simulate experimental constraints. The diameter of the hemispherical head cylindrical impactor is 15 mm and it is simplified into an analytical rigid body to reduce the number of elements and improve computation efficiency. The impactor mass is 1kg, and the impact velocity is set to 3.83 m s −1 to simulate the 7.35 J impact energy.

Element types and mesh density
The eight-node three-dimension cohesive elements (COH3D8) with a thickness 0.0075 mm are inserted between plies to predict the delamination. The maximum damage variable of cohesive elements is set to 0.99 to maintain residual bonding strength. The composite laminate adopts the eight-node solid elements with reduced integration (C3D8R). Considering the damage distribution of the composite material observed in the experiment and the cost of numerical computation, the mesh near the impact area is encrypted, and the mesh size is 1 mm×1 mm. In this paper, there are 35,100 elements in the low-velocity impact finite element numerical model.

Material properties
The material properties of the carbon fibers (Tenax HTS40 12K 300) impregnated with Cycom ® 977-2 epoxy resin are shown in table 2 [21,23,24]. The material properties of the adhesive layer are shown in table 3 [5].

Contact algorithm
The general contact algorithm in ABAQUS/Explicit package is taken into consideration in this paper to simulate the contact between impactor and composite laminate and the self-contact inside the composite laminate. The friction coefficient is realized by a penalty function. According to the results of the 25-27, the friction coefficient between 0°and 0°ply is 0.2, and the friction coefficient between 90°and 90°is 0.8. Some studies take the average as the friction coefficient. Since both the impactor and composite laminates have rough surfaces and the setting of the friction coefficient, which contributes in absorbing impact energy of impactor, has a significant influence on the impact results.In order to obtain the optimal friction coefficient under each failure criterion, this paper sets up multiple sets of friction coefficients for research.

Results and discussion
In this section, the numerical results predicted by 2D Hashin criterion, 3D Hashin criterion and Chang-Chang criterion under three impact energies (7.35 J, 11.03 J, 14.7 J) are compared with the experimental results provided by Shi [21] to validate the finite element model. In addition, the influence of different stacking sequences and different friction coefficient between impactor and composite laminate on the low-velocity impact mechanical response and damage modes are analyzed.

Validation of finite element model
In this part, the difference including energy-time curves, force-time curves, force-displacement curves and damage distribution between FEM and published experimental data provided by Shi [21] is discussed in detailed. The positive influence of s 33 on the computation results is studied.

Global mechanical behavior response
The energy-time curves predicted by each failure criterion under different impact energies are shown in figure 3. The energy-time curves predicted by the three failure criteria are in good agreement with the experimental results. In the initial stage of impact, due to the small damage zone and high stiffness of the composite laminates, the kinetic energy can be quickly absorbed by composite laminate, and most of the energy is converted into internal energy, friction energy and other energy. With the kinetic energy of the impactor decrease, the damage area of the composite laminates continues to increase. As a result, the stiffness of the impact zone gradually decreases, some composite materials even lose their bearing capacity. Therefore, the energy absorption rate decreases, which is manifested as a gradual decrease in the slope of the curve in figure 3. The curve no longer rises, indicating that the kinetic energy is completely absorbed by the composite laminate. At the same time, the composite laminate will release the stored internal energy and react to the impactor, which drives the impactor to rebound. When the rebound speed of the punch no longer increases, the energy absorbed by the composite laminate will tend to a stable value, and the curves in figure 3. will almost no longer change.
By comparison of figures 3(a)-(c), it is found that the prediction accuracy of Chang-Chang criterion and 3D Hashin criterion is higher than that of 2D Hashin criterion. The reason is that the Chang-Chang criterion considers the influence of shear stress in the criterion of fiber failure, and 3D Hashin considers the influence of s 33 on the matrix failure. As the impact energy increases, the prediction results of the 2D Hashin criterion deviate significantly from the experimental results, while the prediction of the 3D Hashin criterion is more accurate. The Chang-Chang criterion will achieve good results under low impact energy (such as 7.35 J). This reason is that the  greater the impact energy, the greater the influence of s 33 on the matrix failure criterion, which makes the damage area of the composite material closer to the experimental results. It is worth noting that the energy absorption value predicted by the three criterion is always smaller than the experimental result. This may be related to some simplification measures in the finite element method, such as simplifying the impactor head to a rigid body.  Figure 4 illustrates the comparison of the numerical results and experimental force-time curves under three impact energies levels. As can be seen, the similar tendency between FEM and experimental data provided by Shi [21]. It is found that the experimental and numerical analysis results suggest that the force-time curve always has a certain vibration during the impact process, and the vibration amplitude is the largest near the peak of the force. The reason is that the damage modes and damage area of composite laminate continue to increase during the impact, resulting in a gradual decrease in the overall stiffness. As a result, the force is continuously redistributed, causing vibration in the force-time curve. The damage modes and damage area of composite laminate barely increase during the rebound stage, resulting in a relatively stable mechanical performance, which makes the force-time curve less vibrating in the rebound stage.
By comparison of figures 4(a)-(c), as the impact energy increases, the value of the maximum force also increases. Among them, the difference between the results predicted by 2D Hashin and the experimental results becomes larger, and the overall fluctuation range is larger. It suggests that for low-speed impact, the effects of s 13 and s 33 on fiber damage and matrix damage cannot be ignored. When the impact energy is low, the trend of force-time curve and value of peak force predicted by the Chang-Chang criterion are in good agreement with the experimental results. However, when the impact energy comes to 14.7 J, the force-time curve predicted by Chang-Chang criteria will deviate from the experimental results during the impact phase. In contrast, the forcetime curve predicted by the 3D Hashin criterion always agrees well with the experimental results, and the forcetime curve in the rebound phase is closer to the experimental curve.
Combining the force-time curve of figure 4 with the displacement-time curve, the force-displacement curve of each failure criterion under three impact energy levels can be obtained, as shown in figure 5. As the impact energy increases, the deformation of the composite laminate becomes larger and larger, which is reflected in the larger displacement value in figures 5(a)-(c). Specific to each failure criterion, it can be found that the maximum deformation predicted by 2D Hashin is always greater than the experimental value, and this gap will increase with the increase of impact energy. For Chang-Chang criterion, the displacement-force curve is in good agreement with the experimental results when subjected to lower impact energy. When the impact energy is 14.7 J, the displacement-force curve starts to deviate from the experimental results. The 3D Hashin criterion always maintains a good prediction accuracy under the three impact energies. It suggests that considering s 33 and shear stress in the damage model can improve the accuracy of finite element prediction, and make the finite element model suitable for predicting the mechanical behavior and damage distribution of composite laminates under various impact energies.

Damage modes
Considering that the matrix damage criteria of the three failure criteria are significantly different, this section will illustrate the distribution of matrix damage and delamination predicted by the three criteria at the impact energy of 14.7 J, as shown in figures 6(a)-(g).
From the figures 6(a)-(d), the overall delamination damage area predicted by the three failure criteria is always larger than the experimental results, and the shape of the delamination damage predicted by 3D Hashin is the closest to the x-ray radiograph. On the other hand, the low-velocity impact finite element model established in this paper can provide conservative suggestions for the impact resistance design. Specifically, it is found that delamination occurred at each interface of the composite laminate, including the completely damaged area (red color spectrum) and the partially damaged area (the rest of the color spectrum). The delamination area of each interface predicted by the 2D Hashin criteria and Chang-Chang criteria is larger than that of 3D Hashin. The reason could be that the damage aera of each single layer predicted by 2D Hashin criteria and Chang-Chang criteria is smaller, which causes the force on each interface to increase. It is found that the delamination propagation direction is mainly parallel to the fiber orientations of the composite single layer below the corresponding interface, which is particularly obvious in the last four interfaces. The phenomenon is consistent with the conclusion given by Zhou [5], Li [8] and Liu [10]. The reason for this phenomenon could be that the delamination propagation direction between the first three interfaces is simultaneously influenced by the impact load and the fiber direction. The impact load will gradually decrease during layer-by-layer transmission, resulting in an increasing influence of fiber direction on the delamination propagation direction in the last four interfaces.
The shapes of the matrix tensile damage predicted by the three failure criteria have good consistency. However, the damage area of each single layer predicted by the 3D Hashin criterion is larger than the other two failure criteria. The reason is that the 3D Hashin criterion considers the influence of s 33 on the matrix tensile damage, so the matrix tensile criterion will be equal to 1 earlier, resulting in a larger damage area. The area of the matrix tensile damage of the last four layers is significantly larger than that of the first four layers. This is because the bending changes of the last four layers are large.
The three failure criteria have the biggest difference in the prediction of matrix compression damage. The 3D Hashin criterion has matrix compression damage on each layer, and the damage area gradually decreases from the first layer to the eighth layer. This is due to the fact that the element deletion technology is not used in the finite element analysis in this paper, which results in the first layer of composite laminate always in contact with the impactor. However, under the prediction of the 2D Hashin and Chang-Chang criteria, the area of matrix compression damage on each single-layer is small, and some single-layer boards do not even show matrix compression damage. This may be because these two damage criteria do not consider the influence of s 33 on the compression damage of the matrix. Figure 7 shows the propagation of fiber tensile damage, matrix tensile damage and matrix compressive damage predicted by 3DHashin criterion under impact energy of 14.7 J, respectively. The moment t=2.9 ms means that the kinetic energy of the impactor has been completely absorbed at this time and it starts to rebound. It is found that the matrix tension damage occurs first in composite laminates and concentrates in the first and eighth layer. The reason is that the first layer is in direct contact with the impactor, resulting in the largest impact load, while the eighth layer has the largest bending deformation. Observing the entire impact process, it can be seen that from the first layer to the eighth layer, the fiber tensile damage and matrix compression damage gradually decrease, while the matrix tensile damage has the opposite trend. This shows that the impact load has a greater influence on fiber tensile damage and matrix compression damage, while the matrix tensile damage is mainly affected by the bending load caused by the deformation of the composite laminate.
Based on the performance of the three failure criteria in the prediction of mechanical behavior and damage distribution, it can be concluded that 3D Hashin is suitable for various impact energies, and the Chang-Chang criterion has better prediction accuracy at low impact energy level.    The reason is that the addition of ±45°plies helps to reduce the overall shear stress level of the composite laminate, which decreases the damage area of the composite laminate.

Influence of friction coefficient
Impact energy is mainly converted into friction energy, internal energy of composite laminates and other forms of energy. The greater the friction coefficient between the impactor and the composite laminate, the greater the friction energy, and the lower the internal energy of the composite laminate, which means that the mechanical behavior of the composite laminate during the impact includes the maximum deformation and the maximum force will also change. Therefore, it is necessary to study the influence of the friction coefficient setting on the computation results. In the literature of Liu [10] and Faggiani [30], the friction coefficient is set to 0.3. In the literature of Zhou [5] and Shi [21], the friction coefficient is set to 0.5. In this section, the mechanical behavior changes of composite laminates after impact are compared under four friction coefficients of 0, 0.3, 0.5, and 0.7. Figure 10 shows the energy-time curves, force-time curves and force-displacement curves using the three failure criteria at the impact energy of 14.7 J. It is found that the increase of the friction coefficient causes the mechanical behavior of the composite laminate to change regularly. In figure 10(a), the curve will deviate significantly from the experimental results by the friction coefficient set of 0, and the error of the predicted energy absorption value will also become larger. The increase of friction coefficient will make the finite element results tend to be consistent and reduce the error. It can be seen from figure 10(b) that with the increase of the friction coefficient, the maximum force predicted by the finite element method will be close to the experimental data provided by Shi [21] and the force-time curves in the rebound phase will be closer to the experimental results. However, when the friction coefficient is increased to 0.5, the prediction accuracy will not be significantly improved. Figure 10(c) shows that increasing the coefficient of friction can effectively reduce the maximum deformation of the composite laminate due to impact. The reason for the above phenomenon is that during the impact process, the kinetic energy of the impactor will be converted into friction energy, the internal energy of the composite laminate and other forms of energy. Increasing the friction coefficient will increase the friction energy and reduce the energy absorbed by the composite material, so the overall strength is improved and closer to the experimental results.
Overall, considering the various positive or negative effects of the friction coefficient between the punch and the composite laminate on the finite element results, when the friction coefficient is set to 0.5-0.7, the prediction accuracy is better.