Progressive multi-damage analysis of composite laminate using higher order zig-zag plate theory

In this article, a finite element model of multi-layer composite laminates with delamination based on higher order zig-zag theory is studied for the progressive delamination analysis. The degrees of freedom are simplified by the continuity conditions of shear stress between layers and the free-surface conditions at the bottom and top of laminates. The numbers of degrees of freedom are not reduced with this method while the model can take the initiation of delamination into consideration in return. Static loading analysis is implemented to simulate the performance in delamination resistance with two different plies and under two different boundary conditions. The present model can present good performance in the prediction of delamination phenomenon.


Introduction
Nowadays, composite materials have been applied extensively in many industrial fields such as aerospace, chemistry, automobile, sport apparatus, and medical instruments, attributing their high specific strength and stiffness, good fatigue and corrosion resistance, and property designability. These merits make them occupy a superiority position in competing against conventional materials. The laminated composite is the most common used composite in various composite structures.
Failure or damage mode of composite material is complex compared with homogeneous material, for example, metal or alloy, attributing their multiphase microstructures. The failure mode of composite laminate is usually classified into three types: fiber failure, matrix failure, and delamination. Failure theories in early stage, for example, Tsai-Wu criterion, 1 are usually phenomenological ones which do not distinguish damage mode in detail. Hashin 2 proposed a well-known failure law in which damage mechanisms for different component phases are distinguished. Puck 3 developed Hashin's work, and their model with the angle of fracture plane can predict unidirectional (UD) composite under transverse compression by introducing a function called the stress exposure factor. A lot of works [4][5][6][7][8][9] have been implemented to verify the reasonability and reliability of their theory numerically and experimentally. In numerical simulation, appropriate material properties degradation should be employed to reflect the residual stiffness or strength of damaged composite corresponding to different failure modes or even for the needs of computing stability. YS Reddy and JN Reddy 10 suggested a stiffness coefficient which resulted good agreements with experiments. Hashin criterion was employed by Hwang and colleagues [11][12][13] and Camanho and Matthews 14 to predict failures of fiber, matrix, and delamination while considering stiffness reduction.
Composite are usually used in thin wall structures; plate and shell theory is an efficient tool for mechanical analysis of these structures. The traditional Kirchhoff thin plate model neglect of transverse shear strain (stress) and normal stress along thickness, which actually have important contributions to delamination damage of composite laminates. For improvements, first-order shear deformation theories (FSDT) 15 and higher order shear deformation theories (HOSDT) were proposed considering the transverse shear stress. Normal stress in the direction of thickness was studied by Pagano 16 and Whitney and Sun. 17 M Aydogdu 18 developed a new displacement mode for laminated plates considering transverse shear deformation and got improved static and dynamic results. Furthermore, Reddy et al. 19 developed a layerwise theory for improving accuracy of multi-layer composite structure. In order to aise computing efficiency, the zig-zag theory was proposed by Lekhnitskii, which has demonstrated conspicuous advantage in hierarchical stress and transverse shear analysis for laminated structures. Zig-zag theories have developed various displacement modes and used that to different problems in recent years. Ambartsumian 20 applied parabolic functions to derive transverse shear stress distribution; YB Cho and RC Averill 21 used a discrete first-order zig-zag theory (FZZT) to analyze laminate and sandwich beam. U Icardi 22 deduced an eight-node zig-zag element; Oh and Cho 23 made an attempt to analyze thermo-electricmechanical coupled smart composite plates with higher order zig-zag theory (HZZT). By adding step functions, zig-zag theories can employed in delamination damage analysis of composite laminates, which can be found in works of U Icardi 24 and RMJ Groh and A Tessler. 25 For delamination analysis, virtual crack closure technique (VCCT), cohesive zone modeling (CZM) and extended finite element method (XFEM) are three of the commonest methods to analyze problems about crack propagation and delamination. VCCT was proposed by Irwin. 26 The VCCT associates the nodal force with displacements around the crack tip. By this method, the strain energy release rate can be easily obtained. Davidson et al. 27 established composite laminate model to study the fracture toughness of mode I and mode II crack with VCCT. Kaczmarek et al. 28 used three-dimensional (3D) VCCT to analyze the delamination at the boundary of laminates. CZM is based on the relative displacement between the upper and lower surface of crack and uses fracture toughness to judge failure. CZM was first proposed by Dugdale. 29 This method includes all three modes of delamination and the constitutive relationship is very simple. Many researchers simulated delamination propagation with CZM and developed CZM. Comanho et al. 30 used CZM to do numerical simulation of delamination propagation. Robinson et al. 31 utilized CZM to deal with the problems about crack propagation under fatigue load. Fan et al. 32 proved that CZM performs well in the analysis of delamination propagation in double cantilever beam (DCB) compared with experiments. The main idea of XFEM is to add additional degrees of freedom to nodes to express the information related to cracks. Motamedi and Mohammadi 33 did the dynamic analysis to laminates with XFEM. Ashari and Mohammadi 34 analyzed the problems about intra-laminar and interlaminar crack. Wang and Waisman 35 developed XFEM with CZM to study delamination propagation.
There are still not much delamination analysis which is applied on HZZT. In this article, a method to process the inter-laminar normal stress is proposed. A displacement mode based on HZZT is presented which can provide continuity of transverse stress and discontinuity of displacement at the interface of layers where delamination happens. The numerical investigation for evolutions of multi-damage including delamination in laminated composite plate was carried out.

Displacement formulation of HZZT
A displacement model of HZZT satisfying two basic conditions, that is, free of transverse on two surfaces of plate and continuity of transverse shear stress along thickness, was proposed by M Cho and JS Kim 36 which can describe multi-delamination of composite laminate. The displacement of this theory is formulated as follow ( Figure 1) The plane where z = 0 is defined as the reference plane. u 0 i (x 1 , x 2 ) represents the displacement in the reference plane, c i (x 1 , x 2 ) indicates the rotation about x i axis, j i and f i are the coefficients of the high-order items in the displacement mode which hardly have directly mechanical meanings, H(x) is the Heaviside function, is the jump of rotation of the (k + 1)th layer compared with that of the kth layer at the kth interface while u k i (x 1 , x 2 ) and w k i (x 1 , x 2 ) reflect the sliding and opening displacement jumps similar to S k i (x 1 , x 2 ), and N is the number of layers. With the condition of free transverse stress at two surfaces, s i3 j z = 0 = 0 and s i3 j z = h = 0, the transverse strain can be deduced as follow On the two surfaces where () , i denotes the partial derivative about x i -coordinate. Due to the continuity of transverse stress between adjacent layers, we have Through combining constitutive relation of laminates, S k i can be associated with material properties and displacement parameters. The relationship is given as follows in a tensor expression where a k ij is only related to material properties. The derivation of a k ij is based on the continuity of transverse shear stress at the interface. The exact process can be obtained as Cho's work. 37 By substituting equation (8) and condition of free transverse shear stress on surfaces into the displacement models (1-3), the displacement fields are formulated as The concrete definition of the coefficients can be found in Appendix 1.

Stiffness matrix of HZZT
In a laminate coordinate system, the resultant forces for each ply in composite laminate relate strains of this ply via stiffness matrix where the superscript (k) denotes the kth layer. Substituting the expressions of strains in Appendix 1 into equation (11), the stresses in kth ply are expressed Integrating the stresses above along the thickness of laminated composite plane, the resultant force of composite is obtained as follows 36 where i are stress resultants. The expressions of elements in stiffness matrixes in equations (14) and (15) are listed in Appendix 2.

Finite element scheme of HZZT
In finite element scheme of HZZT, the displacement vector of a node for an element is where i represents the ith node of the element. The components in equation (16) are defined as with Lagrangian and Hermite interpolations where n is the number of the nodes of one element. In this interpolation, N m is the Lagrangian interpolation functions, and P m and H ym are Hermite interpolation functions. Strain vectors are calculated from nodal displacements of element ð18Þ . .
The expressions of sub-matrixes in the strain-nodal displacement matrixes above are found in Appendix 3.

Multi damages models of composite
Equations (20)-(23) are failure criteria for modes of longitudinal tensile failure, longitudinal compressive failure, transverse tensile failure, and transverse compressive failure, respectively. X T , X C , Y T , Y C , S 23 , S 13 , and S 12 are strengths corresponding to different failure modes under unidirectional load states. It needs to notice that for the transverse compression failure mode, the failure is actually caused by shear stress. And the angle of fracture plane is determined by the transverse normal stresses and shear stress as: As the damage occurs, the engineering constants of the material are declined with the coefficient, as given in Table 1.

Inter-laminate damage modes
The delamination failure criteria is applied as equation (25) where Z int T is the interfacial off-plane normal strength, and S int 13 and S int 23 are the two interfacial transverse shear strength. If all integration points around a node at the interface satisfy equation (25), the interface at this node is considered delaminated. In equation (27), transverse shear stresses s 13 and s 23 are calculated directly using equation (4) as where G 13 and G 23 are two transverse shear moduli.
Ordinary zig-zag model cannot calculate off-plane normal stress s 33 which is important contribution to delamination. Here, HZZT model is developed to calculate off-plane normal stress through introducing equivalent springs at element nodes bonding two adjacent plies, as schematically pictured in Figure 2. The spring stiffness and off-plane normal stress are estimated as where S e is the area of the plate element, and h ply is the distance between the middle planes of the two adjacent plies and it equals the thickness of one ply if each ply has the same thickness in composite laminate. In practical finite element method (FEM) frame, the additional freedom w k of all nodes are released in order to calculate s 33 , and spring stiffness matrix is independent to the stiffness of ordinary stiffness. If debonding (judged by equation (25)) does not happen, the continuities of displacements and shear stresses across the interface are assumed still valid. That means the transverse stresses, as well as the intra-laminate stresses, are calculated in the same way as ordinary zig-zag approach with perfect bonding interface. As soon as the debonding takes place, the stiffness of bonding spring will reduce to zero and the model degenerate to ordinary debonding HZZT one, that is, freedoms of u k i , w k i , and w k , i are all released to realize the delamination occurring or propagation.

Implementation in finite element frame
The damage criteria mentioned above have been implemented in HZZT element on platform of user-defined element (UEL) of general element code Abaqus. If the in-plane damage occurs, the relevant terms in stiffness matrix will be degraded. It needs notice that for an perfect interface between adjacent plies, the freedoms relevant with delamination are forced zero If the delamination occurs between layers, the constrains on u k i , w k i , and w k , i will be released. Figure 3 is the flowchart of progressive damages of HZZT under FEM frame.

Double cantilever beam
To validate the proposed HZZT in this part, simulation of delamination propagation in DCB is conducted. The model is constructed according to Heidari-Rarani and Sayedain. 38 The dimension of the specimen is 150 mm 3 25 mm 3 4:2 mm, and a pre-crack with length of 35 mm is set at one end of the specimen and in the middle of the specimen, which is shown as Figure 4. The size of the element is 1 mm 3 1 mm. The left end of the beam is clamped, and in the right end, a mode I fracture load is applied on the pre-crack, as shown in Figure 4. The material properties of the specimen are given in Table 2.

Orthogonal ply composite laminates
Carbon fiber-reinforced epoxy matrix composite laminates are employed to validate the HZZT approach proposed in this article. The elastic constants and strength parameters of single ply of this composite are listed in Tables 3 and 4, respectively. Two ply patterns are used, that is, (0°, 90°, 90°, 0°) and (45°, -45°, -45°, 45°), which are shown in Figure 5, and they are named ply pattern 1 and ply pattern 2, respectively; the thickness of each ply is 0.25 mm. The dimensions of the two plates of composite laminates are both 101 mm 3 101 mm with the thickness of 1 mm, and they are discretized to square HZZT elements with size 1 mm 3 1 mm. The plates are applied two antisymmetric uniform distributing loads on one side, as illustrated in Figure  6; the loads form a torque which will cause interlaminar shear stress in composite laminate. Two boundary conditions are used: one case is that the side opposite to load-applied edge is clamped and the two other sides are free; the other case is that the three sides except the load-applied edge are all clamped, as illustrated in Figure 6. Table 2. Elastic constants of single ply 1.  Loads with 0.05 N/m increased per step are applied to the composite laminates with different ply patterns and different boundary conditions as Figures 5 and 6 show. The definition of interface is demonstrated in Figure 7.

Mode I interfacial debonding DCB
In this section, load-open displacement curves obtained from DCB model built in section ''Double cantilever beam'' is illustrated in Figure 8, and results from experiment 39 and from other models 38 are also presented in Figure 8 for comparison. The load-open displacement curves show linear pattern during the initial loading stage because the interfacial damage is not severe and the interfacial crack length is considered nearly invariable in this stage. When the delamination propagates, the load declines as the open displacement increases.
As Figure 8 shows, all numerical models fit well with the experimental result in initial stage, but they all over estimate the maximum load capacity compared with experiment measurement. This is explained that for real scenario, interfacial damage occurs earlier than model predictions, and besides, the bridge effects at debonding crack tip might retard the ultimate interface separating. The numerical models, however, do not consider the bridge effect and the complete debonding occurs in sudden way in one element scope. Therefore, the experimental curve exhibits obvious slope reduction before it reaches load peak and it drops from peak point in a more gentle way than model predictions. During the softening stage, the proposed model, and methods of Table 3. Elastic constants of single ply 2.  Table 4. Strength parameters of single ply.  CZM 2D 38 and XFEM-CZM, 38 agrees well with experiment data. Compared with other numerical methods, the HZZT model developed in this study either gives prediction with the same accuracy level of existed models, or it can be built more conveniently.

Multi-mode progressive damages of composite laminates
For case of one-side clamped boundary condition, delamination damage evolutions are illustrated in Figures  9-11. Figure 9 shows the delamination of interface 1, which is counted from bottom surface to top surface as Figure 7 shows. Interfacial damage occurs first at the clamped edge as shown in Figure 9(a), and then it appears at the bending load edge (Figure 9(b)), with load increasing as the existed delamination areas extend and new ones are found in inner region (Figure 9(c)), when load reaches q = 0.9 N/m and almost 85% area of interface 1 fails. Figure 10 is the interfacial damage in the same interface for ply pattern 2. The damage also begins at clamped edge but it extends slower than ply pattern 1. Figure 10(b) shows obviously better performance in delamination resistance than that in ply pattern 1 which is shown in Figure 9 even under much larger load. For the interface 3, however, the delamination at load side is severer than clamped edge as shown in Figure 11(a), then the delamination areas will connect each other with load increasing. For the three sides clamped boundary condition, the interfacial failure of interface 1 for ply pattern 1 starts at the corner of the loading side and then appears at middle of this side, as shown in Figure 12(a) and (b). When load increases, these two delaminated areas extend, connect together with other and isolated interfacial crack areas emerge successively (Figure 12(c) and (d)). Figure 13 shows the delamination evolution in this interface for the ply pattern 2. The delamination is also found at corner of load side (Figures 13(a) and (b)), and subsequently interfacial damage merges at other place of load side ( Figure  13(c)), and then takes place in clamped sides and other isolated regions (Figure 13(d)).

Conclusion
In this article, HZZT is employed and developed for analysis of multi damages (especially for delamination) in composite laminates. It concludes that 1. HZZT is developed to calculate off-plane normal stress by introducing springs at element nodes which bonds with the adjacent plies. The off-plane normal stress is important (for some loading scenario it is the must) for debonding estimation of composite laminate, but however, it is absent in traditional zig-zag models. Therefore, the development made in this study widely extends the application cases of HZZT. 2. The HZZT developed in this study is used to predict load-displacement behavior of a composite DCB, and its prediction agrees well with experimental data and results from other numerical models. Besides, it has advantage in    convenience of model building compared with other existed numerical models. 3. When one side of the laminate is clamped and the opposite is loaded with anti-symmetric uniform load, the delamination usually occurs at the corner connecting clamped border and load free edge initially and then takes place at the loading area. 4. Laminates with three sides clamped usually delaminate at the free side with loads. 5. Delamination areas can spread and intend to join together with the increase of the damage area. 6. Ply pattern in (45°, -45°, -45°, 45°) may have better delamination resistance than that in (0°, 90°, 90°, 0°).

Declaration of conflicting interests
The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding
The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: