Simple Progression Law in Predicting the Damage Onset and Propagation in Composite Notched Laminates

The aim of this article is to simulate the damage initiation and progression in unidirectional (UD) laminates. A three-dimensional (3D) failure criteria of Puck incorporated with degradation scheme is developed. Two types of degradation law known as sudden degradation are used to predict the damage progression in UD laminates. The establishment of constitutive law in progressive damage model (PDM) is achieved through implementation of user subroutines in Abaqus. The failure analysis is applied to various composite stacking sequences and geometries, as well as different fiber reinforced polymer (FRP) composite materials. The comparative studies revealed that the predicted ultimate failure load agree well with those available in the literature. 


Introduction
Damage initiation and evolution in composite structures is very crucial to be analyzed which can occur in several failure mechanisms such as fibre breakage, fibre buckling, matrix cracking, material debonding and also delamination. The failure can occur either individually or several failure modes are combined. Due to increasingly computational capability to handle such complicated behaviors, many attempts are made to predict the failure of composite laminates. Although research concerning failure theories and prediction of strength of composite laminates on fibrereinforced polymer composite has existed for many decades still, no such a robust and reliable method can precisely predict the failure characteristics of the laminates [1]. In general, not all cases fit the experimental results. Failure theories are initially developed for unidirectional composite materials and tailored based on the static regime. Two categories of methods are distinguished: non-interactive and interactive criteria. Non-interactive criteria assume all failure modes are not connected to each other. Failure will occur if stresses/strains in the principal material orientation exceed the respective strengths/strains. Maximum stress/strain criteria imply that technique to identify the maximum failure load. Apart from that, interactive criteria assume interaction between failure mechanisms which focusing on failure surface or envelope. Among the famous interaction criteria that are still used until now are Tsai-Azzi [2], Tsai-Wu [3] and Hoffmann [4]. Based on the above criterion, the laminates are treated as orthotropic materials [5]. Only a single equation is required to relate the interaction between different stress components in the material frame. Hence, still many scientific works adopt it due to their simple forms and acceptability of the accuracy of the strength prediction. A major drawback of those theories is that the approach used not reflects the level of complexity in composite structures. Therefore, the classical failure theories fail to distinguish the mode of failure in composite laminates. Failure mode based theories are needed to overcome that situation and also will make the analysis more meaningful. Hashin and Rotem [6] proposed a two-dimensional (2D) interaction failure initiation theory based on fibre and matrix failure modes, later extended to 3D forms by Hashin [7] itself. Yamada and Sun [8] proposed a prediction method of composite laminate using in-situ shear strength evaluated in the form of cross-ply laminate that suited for the fibre controlled laminates. Motivated by initial work of Yamada and Sun, Chang and Chang [9] developed an extension model by introducing shear non-linearity failure mode. Puck and Schürmann [10] improved the initiation failure theory of Hashin by introducing an angle of fracture plane. Determination of angle of fracture on action plane is governed by shear stresses acting on it. In the case of plane stress, three different modes are distinguished for inter-fiber failure (IFF) modes. Further development of Puck's theories is carried out by Davila et al. [11] which included evaluation of fracture planes for matrix cracking and utilization of fracture mechanics approach to determine the in-situ strength. This criterion is meant for 2D failure analysis and later extended to 3D problems including shear non-linearity developed by Pinho et al. [12] . Another approach to predict the initial failure of composite structures is to use the continuum damage mechanics (CDM) method at ply level. This method is originally developed by Kachanov [13] and is widely used in analyzing the failure of composite laminates. This approach has great ability to predict fibre fracture, matrix cracking and delamination failure mechanism. Ladeveze and Dantec [14] developed an in-plane failure theory for elementary ply using damage mechanics approach. This theory is designed to predict fiber/matrix debonding and matrix micro-cracking in UD composite. Vaziri et al. [15] implemented CDM failure theory on CFRP composite laminates under plane stress condition. They adopted constitutive equations proposed by Matzenmiller et al. [16]. Results showed dependency on mesh size as well as the damage growth on the loading rate. Analysis on ply level is not sufficient to adequately represent the failure of composite structures. Since damage is accumulated before the structure collapse, a proper progression law or post-failure analysis is required to capture the structural response after failure is initiated. The simplest way is to use ply-discount method to reduce the stiffness properties of the material. Once the damage is initiated, the corresponding elastic and strength properties with respect to failure modes are reduced to certain factors or percentages. The new stiffness matrix is determined by using those reduced material properties. The disadvantages of using this technique are mesh dependency and arbitrary degradation parameters. To alleviate the drawbacks mentioned previously, CDM is used by many researchers to imitate the damage evolution of composite structures. Camanho et al. [17] implemented CDM model to control the progression of damage after initiation of failure. Internal variables are related with damage variables depending on the mode of failures. The issue on mesh-dependency is reduced by adopting the crack band model. Also, Barbero et al. [18] developed a damage model suited for individual laminate and later extended its application to laminates. The developed model embedded the theory of CDM and classical thermodynamics. Concept of effective stress is introduced in constitutive equation. Although this degradation technique can reduce the convergence problem during numerical analysis, the difficulty in obtaining parameters determined experimentally became the vital consideration. Failure criteria which can predict the damage initiation and ultimate failure, as well as damage progression are needed to represent the total failure behavior of laminated composite. In general, it is called progressive damage model (PDM). PDM is applied in many structural cases associated to composite laminates such as composite coupon (un-notched), open-hole tension (OHT), openhole compression (OHC), bolted and riveted joints, and many more. The aim of this paper is to establish the damage evaluation procedure of UD composite material using finite element model by emphasizing the following approaches:  Integrate the prominent failure theory of Puck and plydiscounting method as an evaluation tool for composite structures. The selection of this failure criteria is based on feedback obtained from world-wide failure exercise (WWFE)  Predict the onset of failure as well as damage progression using 3D progressive damage model described previously in case of notched laminate composite plate.  Adopt two degradation methods which are used to degrade the stiffness matrix in the occurrence of damage. These approaches play an important role in determining the quality of predicted failure load.  Validate the numerical results with test data reported in literature. The evaluations are carried out by varying the stacking sequences, composite material and diameter of the hole. The significance of degradation phenomena using Abaqus/standard is emphasized.

3D-Damage Model
The ability to anticipate the initiation and progression of damage in composite laminates is crucial to assess their performance properly. Most damage models have at least four failure modes as described by Hashin [7] which are fibre tension, fibre compression, matrix tension and matrix compression. More reliable criteria is required to investigate the failure characteristics of composite parts especially in predicting matrix failure. Thus, Puck's failure theory which is originated from work of Puck and Schürmann [19] is applied in the current analyses based on the excellent performances in estimating the laminated composite failures [5,20,21].
At early development of this theory, maximum stress/strain criteria are used to identify the fibre failure (FF). Later, Puck realized that those criterion are not really physically correct [21], and new analytical equations for unidirectional (UD) composite laminate are written in the following equations: 11  1  12  12  2  3  11   11  1  12  12  2  3  11   1 ...... 0 Where X t and X c are the tensile and compressive strengths of a UD layer in the longitudinal direction; v 12 and v 12f are the major Poisson's ratio for UD lamina and fibre, respectively. The mean stress magnification factor, m σf is assumed to be 1.3 for glass fibre and 1.1 for carbon fibre [23].
For the inter-fibre failure, the concept of fracture plane is introduced since the stresses acting on this plane is assumed to create the fracture. The fracture plane is predicted with respect to the material plane between the angle θ of −90 o and +90 o (180 o due to symmetry plane). General tensor transformations are used in order to obtain the normal and shear stresses acting on the action-plane.
In this analysis, σ n (θ) is the stress normal to fracture plane, τ nt (θ) and τ nl (θ) are the shear stresses parallel and perpendicular to fibre direction in fracture plane, whereby θ is the inclination angle (arbitrary) on fracture plane (see Figure 1). The equation for tensile and compressive inter fibre failure (IFF) are written as: The parameter ψ denotes the shear angle in action plane, R⊥ is failure resistance normal to fibers direction, and R⊥ ψ , R⊥⊥ and R⊥ ‖ are the fracture resistances of the action plane due to the shear stressing. The additional information on Puck's parameters used here is shown in Table 1. Other related equation can be referred in publications [22,20,23].

Damage Evolution
The use of failure criteria independently is not adequate to represent total failure behavior of laminated composite structures. The suitable degradation technique is required to degrade properly the material properties or stiffness matrix itself. The most famous and easy way is to apply ply-discount method whenever damage is identified. This approach is implemented using subroutines and realized using FORTRAN language.

Reducing Stiffness Matrix
To describe the elastic-brittle behavior of fibre-reinforced composites, a constitutive model suited for composite material is used, and a successive attempt is made by Lee et. al [5], and later this approach is called method 1. A 3D-damaged stiffness matrix is written as:   Where C ij is undamaged stiffness component, and G 12 , G 13 and G 23 are the in-plane and out-of-plane shear modulus of composite material. The multiplication factors k 1 , k 2 and k 3 are defined as following: Local damage variables are evaluated by using constant factor as shown in following equation: The subscript i represents ft, fc, mt, and mc. The constant factor successfully implemented by several researchers [5, 24, 25] be-cause of its simplicity of working principle. This approach also known as sudden degradation. The approach is implemented by User-material (UMAT) subroutine. Finally, at the damaged material points, the constitutive model expressed in terms of stressstrain relation can be updated as:

Reducing Elastic Constants
The other approach for degrading the stiffness of material is to reduce the elastic constants accordingly based on failure modes. The progressive model using 3D Puck's formulation is written using user-defined field (USDFLD) subroutine, and is achieved by reducing the elastic properties. This concept is stated as method 2 in further sections. Reduced elastic constants as mentioned below: When fibre fails in tension or compression, the reduced elastic constants are: 11 And for matrix failure in both loading directions, 22 If more than one failure modes occur, all elastic constants are reduced to 1% of its original magnitude, except for Poisson's ratio which are truncated to zero value. The similar pattern of reduction [26] is successfully carried out using different reduction factor. The stresses from previous increment are called into USDFLD subroutine at the beginning of new increment and used to evaluate Puck's criteria. Once criteria is satisfied, the elastic properties are degraded by multiplying the degradation factor with original stiffness values.

Numerical Analyses
To demonstrate the effectiveness of the proposed damage model, numerical analyses are carried out on open-hole coupons of laminated composite subjected to in-plane tensile loading.

Finite Element Model
Finite element (FE) model developed as shown in Figure 3 are performed using Abaqus 6-13, while Figure 2 is used as general dimension for all cases performed in this publication. Even though some of stacking sequences are symmetry, full scale of laminated plies are modelled in order to increase the accuracy of results. The models are simply supported on one side and loaded in longitudinal direction by means of prescribed displacement. Each plies are modelled using reduced integration brick element C3D8R. The plies are stacked using equivalent single layer (ESL) technique provided by Abaqus via lay-up editor. Thus, only one element through thickness is adequate for the simulation of laminated composite plate. To ease the process of extracting loaddisplacement data, a reference point (RF) is located at the end of the composite plate, which the prescribed displacement is applied. Nodes on free face (pulled face) are tied to the RF using equation constrains so that equal displacement of master node (RF) and slave nodes (pulled face) can be achieved. The global element size is 1 mm, whereas smaller element size is modelled at the vicinity of the hole.

User-Subroutines
The requirement of using more sophisticated failure theory is crucially needed to perform failure analysis of composite laminates.
Furthermore, the increase of computational capability has encourage many researchers to implement complicated theory in order to make better prediction. Thus, in this research, a 3D Puck's failure criteria coupled with damage evolution formulation is implemented and realised by the help of user-defined field (USDFLD) and user-material (UMAT) subroutines. The routines are coded using FORTRAN and linked to Abaqus via FORTRAN compiler. The written codes are called for each integration point during analysis of damage in composite structures. Two main contents of the scripts are failure criteria and damage evolution (namely as progressive damage model, PDM). In general, the PDM involved with three main steps: stress/strain evaluation, checking failure initiation, and degrading the elastic constants or stiffness matrix if the damage occurred. The overall process of PDM carried out in this research is displayed in Figure  4.

Computational Effort
The analysis is performed under windows 64-bit platform, Intel Core i7 CPU and 8 GB of random access memory (RAM). The pre-and post-processing are conducted using Abaqus version 6-13, while subroutines are linked using Intel FORTRAN compiler. There were no complicated algorithms which could affect the elapsed time for simple ply discount evolution law. The tangent stiffness operator or DDSDDE (notation in Abaqus subroutine) was equal to damaged stiffness matrix.

Mesh Dependency Analysis
To investigate the influence of mesh size on failure load of each test configurations, a simple composite plate with central hole under tension loading is developed for both method 1 and method 2. The geometry and mesh are illustrated in Figure 5 by modelling size of geometry since the calculation time was very fast. This simple laminate composed of a single layer 90 o (perpendicular to loading direction) and thickness of 1mm. The material properties used was T300/1034-C as described in Table 2 in the next section. The results of analysis at this stage is used to make proper decision on the selection of sufficient number of element needed around the hole. The relationship of mesh dependency is shown in Figure 6 pointed out the change of magnitude of failure load while increasing the number of element (NOE). It can be seen that increasing the number of element from 2688 to 3360, 5040, and 6720 increased the failure load less than 0.6 % for both cases. Thus, the effect of changing the NOE is not significant from 2688 onwards, however showed significant changes less than 2688 elements. From this analysis, it is clearly shown that sufficient NOE was 2688, and are used in all cases for further analyses.  Table 2 are adopted from Lee et al. [5] and Zhao and Zhang [27]. From the evaluation of micromechanical formulation of composite, the out-of-plane shear modulus can be calculated using relation of G 23 = E 22 /2(1+v 23 ). The Christensen's equation of v 23 = v 12 (1-v 12 (E 22 / E 11 ))/(1-v 12 ) is used to calculate the Poisson's ratio in direction 2-3. The experimental data for OHT specimens of T300/1034-C are taken from Lee at al. [5] and Zhang and Zhao [27] , which are originated from work of Chang and Chang [9].  Table 3 shows the comparison of strength between experiments and simulations for two different stacking sequences as well as thickness of T300/1034C composite rectangular plate. Both method 1 and 2 are discussed in terms of accuracy of predicting the ultimate failure load. The most accurate estimation is achieved via the laminate C for method 1 and 2 that give only -3.245% and -5,147% differences, respectively. In general, the approaches implemented here work quite well with test data for error differences ranging from -11.19 % to 6.182 %.
For the sake of comparison, the simulated load-displacement curve for specimens C and D using both methods are shown in Figure 7 in which the quality of predictions are not so close. This is due to different philosophy used for both methods as described in section 2. The effect of giving discount to the stiffness matrix yielded to the premature failure load and displacement, while reducing the elastic properties resulted in over-predicted ultimate load in many cases. However, the selection of degradation factors plays an important role to change the failure estimation in method 2. In the same figure, it is clearly observed that decreasing the thickness of laminate resulted in lower ultimate load and slope of the elastic region. The patterns of each internal damage variables in each ply can be obtained from the simulation. For the specimen A and B, the damage patterns in Figure 8 indicated that the damage path establishes in perpendicular direction with respect to loading direction. Combination of matric cracking and fibre breakage appears in a 0 layer, while matrix failure prominently contributed to failure in a 90 plies. The similar patterns are shown in different sub laminate, thus is not necessary displayed and discussed here. The progression of damage represented fiber and matrix damage for the 0, 45, -45, and 90 layer of the first sub laminate (from top) with thickness of 2.616 mm using method 1 is exhibited in Figure  9, while the other sublaminates are not shown here since the patterns are identical to sublaminate 1. The damage patterns only visualized at the peak load and end of analysis. The failure initiated at the vicinity of hole, in direction perpendicular to loading path. The trends basically followed the nature of damage, where the weakest sections of layers will be the initiation as well as progression of damage path [28]. At 0, 45 and -45 plies, similar damage patterns can be observed for matrix and fibre failure, however, matric cracking is more prominent at 90 layer.    Fig. 9: Failure patterns of specimen B in sublaminate1 using method 1. Table 4 and Figure 10 present the prediction results of failure load and error analysis regarding to three different layup of IM7/977-3 composite specimen. Both prediction methods performed reasonable well in all cases considering the simplicity of progressive damage model applied, although method 1 showed an error of higher than 10 % (-11.42 %) in layup F. Most predicted results are under prediction when compare to experimental data, except for method 2 in specimen G. Overall, method 2 performed better from method 1 in terms of accuracy of the prediction.

Analysis of AS4/PEEK Open-Hole Laminates
The laminated coupons from the third set of data are fabricated from 16 plies of AS4/PEEK composite material [32]. The laminates stacking code was (0/45/90/-45) 2S , whereas the size of the hole varies from the following diameters: D = 2, 3, 5, 8 and 10 mm. The notched laminates had a length of 100 mm and width of 20 mm, while the total thickness was 2 mm. The material properties are collected from many resources (see Table 2). The fibre and matrix damage evolutions at the laminate specimen are shown in Figure 12. Again, because of similarity of damage pattern in other three sublaminates, only the diagrams of four plies represented sublaminate 1 (from outer top) are displayed based on failure modes. Note that only matrix and fibre failures in tension are discussed due to the application of tensile load. The domination of matrix failure in ply 45 o , 90 o , and -45 o exhibited that matrix cracked earlier than fibre breakage. It is also can be seen that fibre damage initiated at the edge of 0 o ply, whereas the damage in the 90 o ply appeared much later than that in the 45 o and -45 o plies. Conclusively, the damage area generated from matrix failure was always larger than those in fibre failures regime.  Table 5 analyses the accuracy of proposed models to predict the ultimate failure by comparing with experimental results for AS4/PEEK laminates. Method 1 produced slightly lower value of failure load when compare to test results, but still within acceptable range. The error about ±5% proved that method 2 performed excellently to predict the peak load for this laminate. Besides, the increase size of diameter from 2 to 10 mm gradually changed the ultimate load as can be seen in Figure 13 Depending on model used, the failure load decreased about 1.5 to 5 kN for different diameter applied in the simulations.

Conclusion
In the present study, the initiation and growth of cracks in fibre reinforced polymer (FRP) composite laminate are satisfactorily estimated using finite element simulations that employ 3D Puck's failure criteria together with ply discounting method as the evolution law. In addition, the proposed models are able to predict the effect of specimen sizes including thickness and diameter of the laminate, as well as different stacking sequences. The performance of method 2 are better than method 1 due to smaller reduction of stiffness matrix. Furthermore, the degradation of stiffness is depend on the selection of degradation factor (arbitrary), which in this case following the rules proposed from literature. In future, the effect of delamination need to be considered in order to predict out of plane failure and deformation. Although some discrepancy between test data and simulations, the proposed PDM's could be considered as a practical design tool for detecting the onset of failure and damage progression in FRP composite laminated structures.