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

Abstract

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.

Keywords

Higher order zig-zag plate theory composite laminate progressive damage analysis delamination

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,¹ are usually phenomenological ones which do not distinguish damage mode in detail. Hashin² proposed a well-known failure law in which damage mechanisms for different component phases are distinguished. Puck³ 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–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¹⁰ suggested a stiffness coefficient which resulted good agreements with experiments. Hashin criterion was employed by Hwang and colleagues^11–13 and Camanho and Matthews¹⁴ 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)¹⁵ 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¹⁶ and Whitney and Sun.¹⁷ M Aydogdu¹⁸ developed a new displacement mode for laminated plates considering transverse shear deformation and got improved static and dynamic results. Furthermore, Reddy et al.¹⁹ 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²⁰ applied parabolic functions to derive transverse shear stress distribution; YB Cho and RC Averill²¹ used a discrete first-order zig-zag theory (FZZT) to analyze laminate and sandwich beam. U Icardi²² deduced an eight-node zig-zag element; Oh and Cho²³ made an attempt to analyze thermo-electric-mechanical 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²⁴ and RMJ Groh and A Tessler.²⁵

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.²⁶ 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.²⁷ established composite laminate model to study the fracture toughness of mode I and mode II crack with VCCT. Kaczmarek et al.²⁸ 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.²⁹ 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.³⁰ used CZM to do numerical simulation of delamination propagation. Robinson et al.³¹ utilized CZM to deal with the problems about crack propagation under fatigue load. Fan et al.³² 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³³ did the dynamic analysis to laminates with XFEM. Ashari and Mohammadi³⁴ analyzed the problems about intra-laminar and inter-laminar crack. Wang and Waisman³⁵ 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.

Theoretical foundation of HZZT

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³⁶ which can describe multi-delamination of composite laminate. The displacement of this theory is formulated as follow (Figure 1)

\begin{matrix} u_{1} (x_{1}, x_{2}, z) = u_{1}^{0} (x_{1}, x_{2}) + ψ_{1} (x_{1}, x_{2}) z \\ + ξ_{1} (x_{1}, x_{2}) z^{2} + ϕ_{1} (x_{1}, x_{2}) z^{3} \\ + \sum_{k = 1}^{N - 1} S_{1}^{k} (x_{1}, x_{2}) (z - z_{k}) H (z - z_{k}) \\ + \sum_{k = 1}^{N - 1} {\bar{u}}_{1}^{k} (x_{1}, x_{2}) H (z - z_{k}) \end{matrix}

(1)

\begin{matrix} u_{2} (x_{1}, x_{2}, z) = u_{2}^{0} (x_{1}, x_{2}) + ψ_{2} (x_{1}, x_{2}) z \\ + ξ_{2} (x_{1}, x_{2}) z^{2} + ϕ_{2} (x_{1}, x_{2}) z^{3} \\ + \sum_{k = 1}^{N - 1} S_{2}^{k} (x_{1}, x_{2}) (z - z_{k}) H (z - z_{k}) \\ + \sum_{k = 1}^{N - 1} {\bar{u}}_{2}^{k} (x_{1}, x_{2}) H (z - z_{k}) \end{matrix}

(2)

u_{3} (x_{1}, x_{2}, z) = w (x_{1}, x_{2}) + \sum_{k = 1}^{N - 1} {\bar{w}}^{k} (x_{1}, x_{2}) H (z - z_{k})

(3)

Figure 1.

Schematic of HZZT: (a) geometry of the model and (b) deformed configurations.

The plane where $z = 0$ is defined as the reference plane. $u_{i}^{0} (x_{1}, x_{2})$ represents the displacement in the reference plane, $ψ_{i} (x_{1}, x_{2})$ indicates the rotation about $x_{i}$ axis, $ξ_{i}$ and $ϕ_{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, $S_{i}^{k} (x_{1}, x_{2})$ is the jump of rotation of the $(k + 1) th$ layer compared with that of the kth layer at the kth interface while ${\bar{u}}_{i}^{k} (x_{1}, x_{2})$ and ${\bar{w}}_{i}^{k} (x_{1}, x_{2})$ reflect the sliding and opening displacement jumps similar to $S_{i}^{k} (x_{1}, x_{2})$ , and $N$ is the number of layers. With the condition of free transverse stress at two surfaces, $σ_{i 3} |_{z = 0} = 0$ and $σ_{i 3} |_{z = h} = 0$ , the transverse strain can be deduced as follow

\begin{matrix} γ_{i 3} = \frac{\partial u_{i}}{\partial z} + \frac{\partial u_{3}}{\partial x_{i}} = ψ_{i} + w_{, i} + 2 ξ_{i} z + 3 ϕ_{i} z^{2} \\ + \sum_{k = 1}^{N - 1} (S_{i}^{k} + {\bar{w}}_{, i}^{k}) H (z - z_{k}) \end{matrix}

(4)

On the two surfaces

γ_{i 3} |_{z = 0} = ψ_{i} + w_{, i} = 0

(5)

\begin{matrix} γ_{i 3} |_{z = h} = ψ_{i} + w_{, i} + 2 ξ_{i} h + 3 ϕ_{i} h^{2} \\ + \sum_{k = 1}^{N - 1} (S_{i}^{k} + {\bar{w}}_{, i}^{k}) = 0 \end{matrix}

(6)

where $()_{, i}$ denotes the partial derivative about x_i-coordinate. Due to the continuity of transverse stress between adjacent layers, we have

τ_{i 3} |_{z = z_{k}^{-}} = τ_{i 3} |_{z = z_{k}^{+}}

(7)

Through combining constitutive relation of laminates, $S_{i}^{k}$ can be associated with material properties and displacement parameters. The relationship is given as follows in a tensor expression

S_{i}^{k} = a_{ij}^{k} ϕ_{j} - {\bar{w}}_{, i}^{k}

(8)

where $a_{ij}^{k}$ is only related to material properties. The derivation of $a_{ij}^{k}$ is based on the continuity of transverse shear stress at the interface. The exact process can be obtained as Cho’s work.³⁷ 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

\begin{matrix} u_{i} = U_{i}^{(0)} + U_{i}^{(1)} z + U_{i}^{(2)} z^{2} + U_{i}^{(3)} z^{3} \\ + \sum_{k = 1}^{N - 1} U_{i}^{k} (z - z_{k}) H (z - z_{k}) + \sum_{k = 1}^{N - 1} {\bar{U}}_{i}^{k} H (z - z_{k}) \end{matrix}

(9)

u_{3} = w + \sum_{k = 1}^{N - 1} {\bar{w}}^{k} H (z - z_{k})

(10)

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

{\begin{matrix} σ_{11}^{(k)} \\ σ_{22}^{(k)} \\ σ_{23}^{(k)} \\ σ_{13}^{(k)} \\ σ_{12}^{(k)} \end{matrix}} = [\begin{matrix} {\bar{Q}}_{1111}^{(k)} & {\bar{Q}}_{1222}^{(k)} & 0 & 0 & {\bar{Q}}_{1612}^{(k)} \\ {\bar{Q}}_{2211}^{(k)} & {\bar{Q}}_{2222}^{(k)} & 0 & 0 & {\bar{Q}}_{2212}^{(k)} \\ 0 & 0 & {\bar{Q}}_{2323}^{(k)} & {\bar{Q}}_{2313}^{(k)} & 0 \\ 0 & 0 & {\bar{Q}}_{1323}^{(k)} & {\bar{Q}}_{1313}^{(k)} & 0 \\ {\bar{Q}}_{1211}^{(k)} & {\bar{Q}}_{1222}^{(k)} & 0 & 0 & {\bar{Q}}_{1212}^{(k)} \end{matrix}] {\begin{matrix} ε_{11}^{(k)} \\ ε_{22}^{(k)} \\ ε_{23}^{(k)} \\ ε_{13}^{(k)} \\ ε_{12}^{(k)} \end{matrix}}

(11)

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

\begin{matrix} σ_{ij}^{(k)} = {\bar{Q}}_{ijmn}^{(k)} \\ [1, z, z^{2}, z^{3}, (z - z_{p}) H (z - z_{p}), H (z - z_{p})] {\begin{matrix} ε_{mn}^{(0)} \\ ε_{mn}^{(1)} \\ ε_{mn}^{(2)} \\ ε_{mn}^{(3)} \\ ε_{mn}^{p} \\ {\bar{ε}}_{mn}^{p} \end{matrix}} \end{matrix}

(12)

σ_{i 3}^{(k)} = {\bar{Q}}_{i 3 m 3}^{(k)} [z, z^{2}, H (z - z_{p})] {\begin{matrix} γ_{m 3}^{(1)} \\ γ_{m 3}^{(2)} \\ γ_{m 3}^{p} \end{matrix}}

(13)

Integrating the stresses above along the thickness of laminated composite plane, the resultant force of composite is obtained as follows³⁶

\begin{matrix} {\begin{matrix} N_{ij} \\ M_{ij} \\ R_{ij}^{(2)} \\ \begin{matrix} R_{ij}^{(3)} \\ {\bar{N}}_{ij}^{l} \\ {\bar{M}}_{ij}^{l} \end{matrix} \end{matrix}} = [\begin{matrix} A_{ijmn}^{(0)} & A_{ijmn}^{(1)} & A_{ijmn}^{(2)} & A_{ijmn}^{(3)} & B_{ijmn}^{k (0)} & E_{ijmn}^{k (0)} \\ A_{ijmn}^{(1)} & A_{ijmn}^{(2)} & A_{ijmn}^{(3)} & A_{ijmn}^{(4)} & B_{ijmn}^{k (1)} & E_{ijmn}^{k (1)} \\ A_{ijmn}^{(2)} & A_{ijmn}^{(3)} & A_{ijmn}^{(4)} & A_{ijmn}^{(5)} & B_{ijmn}^{k (2)} & E_{ijmn}^{k (2)} \\ A_{ijmn}^{(3)} & A_{ijmn}^{(4)} & A_{ijmn}^{(5)} & A_{ijmn}^{(6)} & B_{ijmn}^{k (3)} & E_{ijmn}^{k (3)} \\ B_{ijmn}^{l (0)} & B_{ijmn}^{l (1)} & B_{ijmn}^{l (2)} & B_{ijmn}^{l (3)} & D_{ijmn}^{lk} & F_{ijmn}^{lk} \\ E_{ijmn}^{l (0)} & E_{ijmn}^{l (1)} & E_{ijmn}^{l (2)} & E_{ijmn}^{l (3)} & F_{ijmn}^{lk} & E_{ijmn}^{lk} \end{matrix}] \\ {\begin{matrix} ε_{mn}^{(0)} \\ ε_{mn}^{(1)} \\ ε_{mn}^{(2)} \\ ε_{mn}^{(3)} \\ ε_{mn}^{k} \\ {\bar{ε}}_{mn}^{k} \end{matrix}} \end{matrix}

(14)

{\begin{matrix} V_{i}^{(1)} \\ V_{i}^{(2)} \\ Q_{i}^{l} \end{matrix}} = [\begin{matrix} A_{i 3 m 3}^{(2)} & A_{i 3 m 3}^{(3)} & E_{i 3 m 3}^{k (1)} \\ A_{i 3 m 3}^{(3)} & A_{i 3 m 3}^{(4)} & E_{i 3 m 3}^{k (2)} \\ E_{i 3 m 3}^{l (1)} & E_{i 3 m 3}^{l (2)} & E_{i 3 m 3}^{lk} \end{matrix}] {\begin{matrix} γ_{m 3}^{(1)} \\ γ_{m 3}^{(2)} \\ γ_{m 3}^{k} \end{matrix}}

(15)

where $N_{ij}$ , $M_{ij}$ , $R_{ij}^{(2)}$ , $R_{ij}^{(3)}$ , ${\bar{N}}_{ij}^{l}$ , ${\bar{M}}_{ij}^{l}$ , $V_{i}^{(1)}$ , $V_{i}^{(2)}$ , and $Q_{i}^{l}$ 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

{u_{i}} = {u_{1}^{0}, u_{2}^{0}, w, w_{, 1}, w_{, 2}, ϕ_{1}, ϕ_{2}, {\bar{u}}_{1}^{1}, \dots, {\bar{u}}_{1}^{N - 1}, {\bar{u}}_{2}^{1}, \dots, {\bar{u}}_{2}^{N - 1}, {\bar{w}}^{k}, \dots, {\bar{w}}^{N - 1}, {\bar{w}}_{, 1}^{1}, \dots, {\bar{w}}_{, 1}^{N - 1}, {\bar{w}}_{, 2}^{1}, \dots, {\bar{w}}_{, 2}^{N - 1}}^{T}

(16)

where $i$ represents the $i th$ node of the element.

The components in equation (16) are defined as with Lagrangian and Hermite interpolations

\begin{matrix} [u_{i}^{0}, ϕ_{i}, {\bar{u}}_{i}^{k}] = \sum_{m = 1}^{n} N_{m} [{(u_{i}^{0})}_{m}, {(ϕ_{i})}_{m}, {({\bar{u}}_{i}^{k})}_{m}] \\ w = \sum_{m = 1}^{n} [P_{m} {(w)}_{m} + H_{xm} {(w_{, x})}_{m} + H_{ym} {(w_{, y})}_{m}] \\ {\bar{w}}^{i} = \sum_{m = 1}^{n} [P_{m} {({\bar{w}}^{i})}_{m} + H_{xm} {({\bar{w}}_{, x}^{i})}_{m} + H_{ym} {({\bar{w}}_{, y}^{i})}_{m}] \end{matrix}

(17)

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

\begin{matrix} {ε_{ij}} = {[B]}_{b} {u_{n}} = \sum_{m = 1}^{n} [B_{mb, ij}] {u_{m}} \\ = [B_{1 b, ij}, B_{2 b, ij}, \dots, B_{nb, ij}] {\begin{matrix} {u_{1}} \\ {u_{2}} \\ ⋮ \\ {u_{n}} \end{matrix}} \end{matrix}

(18)

\begin{matrix} {γ_{i 3}} = {[B]}_{s} {u_{n}} = \sum_{m = 1}^{n} [B_{ms, ij}] {u_{m}} \\ = [B_{1 s, ij}, B_{2 s, ij}, \dots, B_{ns, ij}] {\begin{matrix} {u_{1}} \\ {u_{2}} \\ ⋮ \\ {u_{n}} \end{matrix}} \end{matrix}

(19)

The expressions of sub-matrixes in the strain-nodal displacement matrixes above are found in Appendix 3.

Multi damages models of composite

Intra-laminate damage modes

3D Hashin criterion is employed to predict the initiation of damage

f_{ft} = {(\frac{σ_{11}}{X_{T}})}^{2} + \frac{σ_{12}^{2} + σ_{13}^{2}}{S_{12}^{2}} \geq 1 (σ_{11} \geq 0)

(20)

f_{fc} = {(\frac{σ_{11}}{X_{C}})}^{2} \geq 1 (σ_{11} < 0)

(21)

\begin{matrix} f_{mt} = {(\frac{σ_{22} + σ_{33}}{Y_{T}})}^{2} + \frac{σ_{23}^{2} - σ_{22} σ_{33}}{S_{23}^{2}} \\ + \frac{σ_{12}^{2} + σ_{13}^{2}}{S_{12}^{2}} \geq 1 (σ_{22} + σ_{33} \geq 0) \end{matrix}

(22)

\begin{matrix} f_{mc} = ({(\frac{Y_{C}}{2 S_{23}})}^{2} - 1) (\frac{σ_{22} + σ_{33}}{Y_{C}}) \\ + {(\frac{σ_{22} + σ_{33}}{2 S_{23}})}^{2} + \frac{σ_{23}^{2} - σ_{22} σ_{33}}{S_{23}^{2}} \\ + \frac{σ_{12}^{2} + σ_{13}^{2}}{S_{12}^{2}} \geq 1 (σ_{22} + σ_{33} < 0) \end{matrix}

(23)

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:

α = \frac{1}{2} \arctan^{- 1} \frac{2 τ_{23}}{σ_{22} - σ_{33}}

(24)

As the damage occurs, the engineering constants of the material are declined with the coefficient, as given in Table 1.

Table 1.

Stiffness degradation under different failure modes.

Mode	$E_{1}$	$E_{2}$	$E_{3}$	$G_{12}$	$G_{23}$	$G_{13}$	$v_{12}$	$v_{23}$	$v_{13}$
$ft$	0.01	–	–	0.01	–	0.01	0.01	–	0.01
$fc$	0.2	–	–	0.01	0.01	0.01	0.01	–	0.01
$mt$	–	0.01	–	0.01	–	0.01	–	–	–
$mc$	–	0.2	–	0.2	0.2	–	0.2	–	–

Inter-laminate damage modes

The delamination failure criteria is applied as equation (25)

f_{df} = {(\frac{σ_{33}}{Z_{T}^{int}})}^{2} + {(\frac{σ_{13}}{S_{13}^{int}})}^{2} + {(\frac{σ_{23}}{S_{23}^{int}})}^{2} \geq 1 (σ_{22} + σ_{33} \geq 0)

(25)

where $Z_{T}^{int}$ is the interfacial off-plane normal strength, and $S_{13}^{int}$ and $S_{23}^{int}$ 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 $σ_{13}$ and $σ_{23}$ are calculated directly using equation (4) as

σ_{13} = G_{13} γ_{13}, σ_{13} = G_{13} γ_{13}

(26)

where $G_{13}$ and $G_{23}$ are two transverse shear moduli. Ordinary zig-zag model cannot calculate off-plane normal stress $σ_{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

K_{s} = E_{2} \frac{S_{e}}{4}

(27)

σ_{33} = E_{2} \frac{{\bar{w}}^{k}}{h_{ply}}

(28)

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 ${\bar{w}}^{k}$ of all nodes are released in order to calculate $σ_{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 ${\bar{u}}_{i}^{k}$ , ${\bar{w}}_{i}^{k}$ , and ${\bar{w}}_{, i}^{k}$ are all released to realize the delamination occurring or propagation.

Figure 2.

Model of off-plane spring bonding adjacent plies.

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 ${\bar{u}}_{i}^{k} = {\bar{w}}_{i}^{k} = {\bar{w}}_{, i}^{k} = 0$ . If the delamination occurs between layers, the constrains on ${\bar{u}}_{i}^{k}$ , ${\bar{w}}_{i}^{k}$ , and ${\bar{w}}_{, i}^{k}$ will be released. Figure 3 is the flowchart of progressive damages of HZZT under FEM frame.

Figure 3.

Flowchart of progressive damages of HZZT.

Numerical models

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

The dimension of the specimen is $150 mm \times 25 mm \times 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 \times 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.

Figure 4.

Dimension and boundary (loading) conditions of the DCB specimen.

Table 2.

Elastic constants of single ply 1.

$E_{1}$ $(GPa)$	$E_{2}$ $(GPa)$	$G_{12}$ $(GPa)$	$v_{12}$	$Z_{T}^{int}$ $(MPa)$	$S_{13}^{int}$ $(MPa)$	$S_{23}^{int}$ $(MPa)$
33.5	10.23	4.26	0.27	15	92	92

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 × 101 mm with the thickness of 1 mm, and they are discretized to square HZZT elements with size 1 mm × 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 3.

Elastic constants of single ply 2.

$E_{1}$ $(GPa)$	$E_{2}$ $(GPa)$	$E_{3}$ $(GPa)$	$G_{12}$ $(GPa)$	$G_{13}$ $(GPa)$	$G_{23}$ $(GPa)$	$v_{12}$	$v_{13}$	$v_{23}$
157.2	8.96	8.96	5.08	5.08	2.99	0.32	0.32	0.5

Table 4.

Strength parameters of single ply.

$X_{T}$ $(MPa)$	$X_{C}$ $(MPa)$	$Y_{T}$ $(MPa)$	$Y_{C}$ $(MPa)$	$Z_{T}$ $(MPa)$	$S_{12}$ $(MPa)$	$S_{13}$ $(MPa)$	$S_{23}$ $(MPa)$	$G_{IC}$ $(J / m^{2})$	$G_{IIC}$ $(J / m^{2})$
2325.6	1199.7	62.0	199.7	62.0	92.4	92.4	92.4	240	739

Figure 5.

Ply patterns of the composite plate: (a) ply pattern 1 and (b) ply pattern 2.

Figure 6.

Loading and boundary conditions of the composite plate: (a) one edge clamped and (b) three edges clamped.

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.

Figure 7.

Definition of interface.

Results and discussion

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³⁹ and from other models³⁸ 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.

Figure 8.

Comparison of load-open displacement curves of DCB specimen from proposed model, experiment and other models.

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 CZM 2D³⁸ and XFEM-CZM,³⁸ 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)).

Figure 9.

Delamination development of interface 1 in one side clamped composite laminate of layer pattern 1: (a) q = 0.35 N/m, (b) q = 0.45 N/m, (c) q = 0.7 N/m, and (d) q = 0.9 N/m.

Figure 10.

Delamination development of interface 1 of one side clamped composite laminate of layer pattern 2: (a) q = 0.45 N/m and (b) q = 1.5 N/m.

Figure 11.

Delamination development of interface 3 in one side clamped composite laminate of layer pattern 2: (a) q = 0.7 N/m and (b) q = 1.5 N/m.

Figure 12.

Delamination development of interface 1 in three sides clamped composite laminate of layer pattern 1: (a) q = 0.35 N/m, (b) q = 0.45 N/m, (c) q = 1.2 N/m, and (d) q = 1.3 N/m.

Figure 13.

Delamination development of interface 1 in three sides clamped composite laminate of layer pattern 2: (a) q = 0.35 N/m, (b) q = 0.45 N/m, (c) q = 0.8 N/m, and (d) q = 1.3 N/m.

Conclusion

In this article, HZZT is employed and developed for analysis of multi damages (especially for delamination) in composite laminates. It concludes that

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.

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.

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.

Laminates with three sides clamped usually delaminate at the free side with loads.

Delamination areas can spread and intend to join together with the increase of the damage area.

Ply pattern in (45°, –45°, –45°, 45°) may have better delamination resistance than that in (0°, 90°, 90°, 0°).

Footnotes

Appendix 1

The coefficients in equations (9) and (10) are

\begin{matrix} U_{i}^{(0)} = u_{i}^{0}, U_{i}^{(1)} = - w_{, i}, U_{i}^{(2)} = - \frac{3 h}{2} ϕ_{i} - \frac{z^{2}}{2 h} \sum_{k = 1}^{N - 1} a_{ij}^{k} ϕ_{j} \\ U_{i}^{(3)} = ϕ_{i}, U_{i}^{k} = a_{ij}^{k} ϕ_{j} - {\bar{w}}_{, i}^{k}, {\bar{U}}_{i}^{k} = {\bar{u}}_{i}^{k} \end{matrix}

And according to constitutive equations, the strain tensor can be shown as

\begin{matrix} ε_{ij} = ε_{ij}^{(0)} + ε_{ij}^{(1)} z + ε_{ij}^{(2)} z^{2} + ε_{ij}^{(3)} z^{3} \\ + \sum_{k = 1}^{N - 1} [ε_{ij}^{k} (z - z_{k}) + {\bar{ε}}_{ij}^{k}] H (z - z_{k}) \\ γ_{i 3} = γ_{i 3}^{(1)} z + γ_{i 3}^{(2)} z^{2} + \sum_{k = 1}^{N - 1} γ_{i 3}^{k} H (z - z_{k}) \end{matrix}

where

\begin{matrix} ε_{ij}^{(0)} = \frac{1}{2} (u_{i, j}^{0} + u_{j, i}^{0}), ε_{ij}^{(1)} = - \frac{1}{2} (w_{, ij} + w_{, ji}) \\ ε_{ij}^{(2)} = - \frac{3 h}{4} (ϕ_{i, j} + ϕ_{j, i}) - \frac{1}{4 h} \sum_{k = 1}^{N - 1} (a_{im}^{k} ϕ_{m, j} + a_{jn}^{k} ϕ_{n, i}) \\ ε_{ij}^{(3)} = \frac{1}{2} (ϕ_{i, j} + ϕ_{j, i}), ε_{ij}^{k} = \frac{1}{2} (a_{im}^{k} ϕ_{m, j} + a_{jn}^{k} ϕ_{n, i} - 2 {\bar{w}}_{, ij}^{k}) \\ {\bar{ε}}_{ij}^{k} = \frac{1}{2} ({\bar{u}}_{i, j}^{k} + {\bar{u}}_{j, i}^{k}), γ_{i 3}^{(1)} = - (3 h ϕ_{i} + \frac{1}{h} \sum_{k = 1}^{N - 1} a_{im}^{k} ϕ_{m}) \\ γ_{i 3}^{(2)} = 3 ϕ_{i}, γ_{i 3}^{k} = a_{im}^{k} ϕ_{m} \end{matrix}

Appendix 2

The expressions of elements in stiffness matrixes in equations (14) and (15) are

\begin{matrix} A_{ijmn}^{(p)} = \sum_{k = 0}^{N - 1} \int_{z_{k}}^{z_{k + 1}} {\bar{Q}}_{ijmn}^{k + 1} z^{p} dz, (p = 0, 1, 2, 3, 4, 5, 6) \\ B_{ijmn}^{l (p)} = \sum_{k = 0}^{N - 1} \int_{z_{k}}^{z_{k + 1}} {\bar{Q}}_{ijmn}^{k + 1} (z - z_{l}) H (z - z_{l}) z^{p} dz, (p = 0, 1, 2, 3) \\ E_{ijmn}^{l (p)} = \sum_{k = 0}^{N - 1} \int_{z_{k}}^{z_{k + 1}} {\bar{Q}}_{ijmn}^{k + 1} H (z - z_{l}) z^{p} dz, (p = 0, 1, 2, 3) \\ E_{ijmn}^{lk} = \sum_{k = 0}^{N - 1} \int_{z_{k}}^{z_{k + 1}} {\bar{Q}}_{ijmn}^{k + 1} H (z - z_{l}) H (z - z_{k}) dz, \\ F_{ijmn}^{lk} = \sum_{k = 0}^{N - 1} \int_{z_{k}}^{z_{k + 1}} {\bar{Q}}_{ijmn}^{k + 1} (z - z_{l}) H (z - z_{l}) H (z - z_{k}) dz, \\ D_{ijmn}^{lk} = \sum_{k = 0}^{N - 1} \int_{z_{k}}^{z_{k + 1}} {\bar{Q}}_{ijmn}^{k + 1} (z - z_{l}) (z - z_{k}) H (z - z_{l}) H (z - z_{k}) dz, \\ A_{i 3 m 3}^{(p)} = \sum_{k = 0}^{N - 1} \int_{z_{k}}^{z_{k + 1}} {\bar{Q}}_{ii 3 m 3}^{k + 1} z^{p} dz, (p = 2, 3, 4) \\ E_{i 3 m 3}^{l (p)} = \sum_{k = 0}^{N - 1} \int_{z_{k}}^{z_{k + 1}} {\bar{Q}}_{i 3 m 3}^{k + 1} H (z - z_{l}) z^{p} dz, (p = 1, 2) \\ E_{i 3 m 3}^{lk} = \sum_{k = 0}^{N - 1} \int_{z_{k}}^{z_{k + 1}} {\bar{Q}}_{i 3 m 3}^{k + 1} H (z - z_{l}) H (z - z_{k}) dz . \end{matrix}

Appendix 3

The expressions of sub-matrixes in the strain-nodal displacement matrixes in equations (19) and (20) are

\begin{matrix} [B_{mb, 11}] = [B_{mb, 11}^{u}, B_{mb, 11}^{d}] \\ [B_{mb, 22}] = [B_{mb, 22}^{u}, B_{mb, 22}^{d}] \\ [B_{mb, 12}] = [B_{mb, 12}^{u}, B_{mb, 12}^{d}] \end{matrix}

[B_{mb, 11}^{u}] = [\begin{matrix} N_{m, x} & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & - P_{m, xx} & - H_{xm, xx} & - H_{ym, xx} & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & - (\frac{3 h}{2} + \frac{1}{2 h} A_{11}) N_{m, x} & - \frac{1}{2 h} A_{12} N_{m, x} \\ 0 & 0 & 0 & 0 & 0 & N_{m, x} & 0 \\ 0 & 0 & 0 & 0 & 0 & a_{11}^{1} N_{m, x} & a_{12}^{1} N_{m, x} \\ ⋮ & ⋮ & ⋮ & ⋮ & ⋮ & ⋮ & ⋮ \\ 0 & 0 & 0 & 0 & 0 & a_{11}^{N - 1} N_{m, x} & a_{12}^{N - 1} N_{m, x} \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 \end{matrix}]

[B_{mb, 11}^{d}] = [\begin{matrix} {[0]}_{1 \times (N - 1)} & {[0]}_{1 \times (N - 1)} & {[0]}_{1 \times (N - 1)} & {[0]}_{1 \times (N - 1)} & {[0]}_{1 \times (N - 1)} \\ {[0]}_{1 \times (N - 1)} & {[0]}_{1 \times (N - 1)} & {[0]}_{1 \times (N - 1)} & {[0]}_{1 \times (N - 1)} & {[0]}_{1 \times (N - 1)} \\ {[0]}_{1 \times (N - 1)} & {[0]}_{1 \times (N - 1)} & {[0]}_{1 \times (N - 1)} & {[0]}_{1 \times (N - 1)} & {[0]}_{1 \times (N - 1)} \\ {[0]}_{1 \times (N - 1)} & {[0]}_{1 \times (N - 1)} & {[0]}_{1 \times (N - 1)} & {[0]}_{1 \times (N - 1)} & {[0]}_{1 \times (N - 1)} \\ {[0]}_{(N - 1) \times (N - 1)} & {[0]}_{(N - 1) \times (N - 1)} & {[- P_{m, xx}]}_{(N - 1) \times (N - 1)} & {[- H_{xm, xx}]}_{(N - 1) \times (N - 1)} & {[- H_{ym, xx}]}_{(N - 1) \times (N - 1)} \\ {[N_{m, x}]}_{(N - 1) \times (N - 1)} & {[0]}_{(N - 1) \times (N - 1)} & {[0]}_{(N - 1) \times (N - 1)} & {[0]}_{(N - 1) \times (N - 1)} & {[0]}_{(N - 1) \times (N - 1)} \end{matrix}]

[B_{mb, 22}^{u}] = [\begin{matrix} 0 & N_{m, y} & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & - P_{m, yy} & - H_{xm, yy} & - H_{ym, yy} & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & - \frac{1}{2 h} A_{21} N_{m, y} & - (\frac{3 h}{2} + \frac{1}{2 h} A_{22}) N_{m, y} \\ 0 & 0 & 0 & 0 & 0 & 0 & N_{m, y} \\ 0 & 0 & 0 & 0 & 0 & a_{21}^{1} N_{m, y} & a_{22}^{1} N_{m, y} \\ ⋮ & ⋮ & ⋮ & ⋮ & ⋮ & ⋮ & ⋮ \\ 0 & 0 & 0 & 0 & 0 & a_{21}^{N - 1} N_{m, y} & a_{22}^{N - 1} N_{m, y} \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 \end{matrix}]

[B_{mb, 22}^{d}] = [\begin{matrix} {[0]}_{1 \times (N - 1)} & {[0]}_{1 \times (N - 1)} & {[0]}_{1 \times (N - 1)} & {[0]}_{1 \times (N - 1)} & {[0]}_{1 \times (N - 1)} \\ {[0]}_{1 \times (N - 1)} & {[0]}_{1 \times (N - 1)} & {[0]}_{1 \times (N - 1)} & {[0]}_{1 \times (N - 1)} & {[0]}_{1 \times (N - 1)} \\ {[0]}_{1 \times (N - 1)} & {[0]}_{1 \times (N - 1)} & {[0]}_{1 \times (N - 1)} & {[0]}_{1 \times (N - 1)} & {[0]}_{1 \times (N - 1)} \\ {[0]}_{1 \times (N - 1)} & {[0]}_{1 \times (N - 1)} & {[0]}_{1 \times (N - 1)} & {[0]}_{1 \times (N - 1)} & {[0]}_{1 \times (N - 1)} \\ {[0]}_{(N - 1) \times (N - 1)} & {[0]}_{(N - 1) \times (N - 1)} & {[- P_{m, yy}]}_{(N - 1) \times (N - 1)} & {[- H_{xm, yy}]}_{(N - 1) \times (N - 1)} & {[- H_{ym, yy}]}_{(N - 1) \times (N - 1)} \\ {[0]}_{(N - 1) \times (N - 1)} & {[N_{m, y}]}_{(N - 1) \times (N - 1)} & {[0]}_{(N - 1) \times (N - 1)} & {[0]}_{(N - 1) \times (N - 1)} & {[0]}_{(N - 1) \times (N - 1)} \end{matrix}]

[B_{mb, 12}^{u}] = [\begin{matrix} \frac{1}{2} N_{m, y} & \frac{1}{2} N_{m, x} & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & - P_{m, xy} & - H_{xm, xy} & - H_{ym, xy} & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & - \frac{1}{4 h} A_{21} N_{m, x} - \frac{3 h}{4} N_{m, y} - \frac{1}{4 h} A_{11} N_{m, y} & - \frac{1}{4 h} A_{12} N_{i, y} - \frac{3 h}{4} N_{i, x} - \frac{1}{4 h} A_{22} N_{i, x} \\ 0 & 0 & 0 & 0 & 0 & \frac{1}{2} N_{m, y} & \frac{1}{2} N_{m, x} \\ 0 & 0 & 0 & 0 & 0 & \frac{1}{2} a_{21}^{1} N_{m, x} + \frac{1}{2} a_{11}^{1} N_{m, y} & \frac{1}{2} a_{12}^{1} N_{m, y} + \frac{1}{2} a_{22}^{1} N_{m, x} \\ ⋮ & ⋮ & ⋮ & ⋮ & ⋮ & ⋮ & ⋮ \\ 0 & 0 & 0 & 0 & 0 & \frac{1}{2} a_{21}^{N - 1} N_{m, x} + \frac{1}{2} a_{11}^{N - 1} N_{m, y} & \frac{1}{2} a_{12}^{N - 1} N_{m, y} + \frac{1}{2} a_{22}^{N - 1} N_{m, x} \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 \end{matrix}]

[B_{mb, 12}^{d}] = [\begin{matrix} {[0]}_{1 \times (N - 1)} & {[0]}_{1 \times (N - 1)} & {[0]}_{1 \times (N - 1)} & {[0]}_{1 \times (N - 1)} & {[0]}_{1 \times (N - 1)} \\ {[0]}_{1 \times (N - 1)} & {[0]}_{1 \times (N - 1)} & {[0]}_{1 \times (N - 1)} & {[0]}_{1 \times (N - 1)} & {[0]}_{1 \times (N - 1)} \\ {[0]}_{1 \times (N - 1)} & {[0]}_{1 \times (N - 1)} & {[0]}_{1 \times (N - 1)} & {[0]}_{1 \times (N - 1)} & {[0]}_{1 \times (N - 1)} \\ {[0]}_{1 \times (N - 1)} & {[0]}_{1 \times (N - 1)} & {[0]}_{1 \times (N - 1)} & {[0]}_{1 \times (N - 1)} & {[0]}_{1 \times (N - 1)} \\ {[0]}_{(N - 1) \times (N - 1)} & {[0]}_{(N - 1) \times (N - 1)} & {[- P_{m, xy}]}_{(N - 1) \times (N - 1)} & {[- H_{xm, xy}]}_{(N - 1) \times (N - 1)} & {[- H_{ym, xy}]}_{(N - 1) \times (N - 1)} \\ {[\frac{1}{2} N_{m, y}]}_{(N - 1) \times (N - 1)} & {[\frac{1}{2} N_{m, x}]}_{(N - 1) \times (N - 1)} & {[0]}_{(N - 1) \times (N - 1)} & {[0]}_{(N - 1) \times (N - 1)} & {[0]}_{(N - 1) \times (N - 1)} \end{matrix}]

[B_{sm, 13}] = [\begin{matrix} 0 & 0 & 0 & 0 & 0 & (- 3 h - \frac{1}{h} A_{11}) N_{m} & - \frac{1}{h} A_{12} N_{m} \\ 0 & 0 & 0 & 0 & 0 & 3 N_{m} & 0 \\ 0 & 0 & 0 & 0 & 0 & a_{11}^{1} N_{m} & a_{12}^{1} N_{m} \\ ⋮ & ⋮ & ⋮ & ⋮ & ⋮ & ⋮ & ⋮ \\ 0 & 0 & 0 & 0 & 0 & a_{11}^{N - 1} N_{m} & a_{12}^{N - 1} N_{m} \end{matrix} {[0]}_{(N + 1) \times 5 (N + 1)}]

[B_{sm, 23}] = [\begin{matrix} 0 & 0 & 0 & 0 & 0 & - \frac{1}{h} A_{21} N_{m} & (- 3 h - \frac{1}{h} A_{22}) N_{m} \\ 0 & 0 & 0 & 0 & 0 & 0 & 3 N_{m} \\ 0 & 0 & 0 & 0 & 0 & a_{21}^{1} N_{m} & a_{22}^{1} N_{m} \\ ⋮ & ⋮ & ⋮ & ⋮ & ⋮ & ⋮ & ⋮ \\ 0 & 0 & 0 & 0 & 0 & a_{21}^{N - 1} N_{m} & a_{22}^{N - 1} N_{m} \end{matrix} {[0]}_{(N + 1) \times 5 (N + 1)}]

and

A_{ij} = \sum_{k = 1}^{N - 1} a_{ij}^{k}

\begin{matrix} {[N_{m, i}]}_{(N - 1) \times (N - 1)} = [\begin{matrix} N_{m, i} \\ ⋱ \\ N_{m, i} \end{matrix}] \\ {[P_{m, ij}]}_{(N - 1) \times (N - 1)} = [\begin{matrix} P_{m, ij} \\ ⋱ \\ P_{m, ij} \end{matrix}] \\ {[H_{lm, ij}]}_{(N - 1) \times (N - 1)} = [\begin{matrix} H_{lm, ij} \\ ⋱ \\ H_{lm, ij} \end{matrix}] \end{matrix}

Handling Editor: José Correia

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: This work was supported by Postgraduate Research & Practice Innovation Program of Jiangsu Province (Grant No. KYCX18_0251), National Natural Science Foundation of China (Grant Nos 11872205, 11272147), SKL Open Fund (Grant No. MCMS-0218G01), and Priority Academic Program Development of Jiangsu Higher Education Institutions.

ORCID iD

Zhaotian Wu

References

Tsai

. A general theory of strength for anisotropic materials. J Compos Mater 1971; 5: 58–80.

Hashin

. Failure criteria for unidirectional fiber composites. J Appl Mech 1980; 47: 329–334.

Puck

. Calculating the strength of glass fiber/plastic laminates under combined load. Ger Plast 1969; 55: 18–19.

Puck

Kopp

Knops

. Guidelines for the determination of the parameters in Puck’s action plane strength criterion. Compos Sci Technol 2002; 62: 371–378.

Pachajoa

Frances

Lee

. Stress and failure analysis of composite structures. Eng Fract Mech 1995; 50: 883–902.

Rakesh

Singh

Kumar

. Failure prediction in glass fiber reinforced plastics laminates with drilled hole under uni-axial loading. Mater Des 2010; 31: 3002–3007.

Liu

Xing

Zheng

. Failure analysis of carbon fiber/epoxy composite cylindrical laminates using explicit finite element method. Composites Part B 2014; 56: 54–61.

Schirmaier

Weiland

Kärger

, et al. A new efficient and reliable algorithm to determine the fracture angle for Puck’s 3D matrix failure criterion for UD composites. Compos Sci Technol 2014; 100: 19–25.

Reinoso

Catalanotti

Blázquez

, et al. A consistent anisotropic damage model for laminated fiber-reinforced composites using the 3D-version of the Puck failure criterion. Int J Solids Struct 2017; 126–127: 37–53.

10.

Reddy

. Three-dimensional finite element progressive failure analysis of composite laminates under axial extension. J Compos Tech Res 1993; 15: 73–87.

11.

Hwang

Sun

. Failure analysis of laminated composites by using iterative three-dimensional finite element method. Compos Struct 1989; 33: 41–47.

12.

Hung

Chang

. Bearing failure of bolted composite joints. J Compos Mater 1996; 30: 1359–1400.

13.

Hung

Chang

. Strength envelope of bolted composite joints under bypass loads. J Compos Mater 1996; 30: 1402–1435.

14.

Camanho

Matthews

. A progressive damage model for mechanically fastened joints in composite laminates. J Compos Mater 1999; 33: 2248–2280.

15.

Reddy

. A generalization of two-dimensional theories of laminated composite plates. Commun Appl Numer Methods 1987; 3: 173–180.

16.

Pagano

. Mechanics of composite materials. Oxford: Pergamon, 1973, pp. 268–286.

17.

Whitney

Sun

. A higher order theory for extensional motion of laminated composites. J Sound Vib 1973; 30: 85–97.

18.

Aydogdu

. A new shear deformation theory for laminated composite plates. Compos Struct 2009; 89: 94–101.

19.

Reddy

Barbero

Teply

. A plate bending element based on a generalized laminate plate theory. Int J Numer Meth Eng 1989; 28: 2275–2292.

20.

Ambartsumian

. On a general theory of anisotropic shells. J Appl Math Mech 1958; 22: 305–319.

21.

Cho

Averill

. An improved theory and finite-element model for laminated composite and sandwich beams using first-order zig-zag sublaminate approximations. Compos Struct 1997; 37: 281–298.

22.

Icardi

. Eight-noded zig-zag element for deflection and stress analysis of plates with general lay-up. Composites Part B 1998; 29: 425–441.

23.

Cho

. Higher order zig-zag theory for fully coupled thermo-electric-mechanical smart composite plates. J Korean Soc Aeronaut Space Sci 2002; 30: 9–14.

24.

Icardi

Zardo

. C₀ plate element for delamination damage analysis, based on a zig-zag model and strain energy updating. Int J Impact Eng 2005; 31: 579–606.

25.

Groh

RMJ

Tessler

. Computationally efficient beam elements for accurate stresses in sandwich laminates and laminated composites with delaminations. Comput Methods Appl Mech Eng 2017; 320: 369–395.

26.

Irwin

. Analysis of stresses and strains near the end of a crack traversing a plate. Appl Mech 1957; 24: 361–364.

27.

Davidson

Kruger

Konig

. Three-dimensional analysis of center-delaminated unidirectional and multidirectional single-leg bending specimens. Compos Sci Technol 1995; 54: 385–394.

28.

Kaczmarek

Wisnom

Jones

. Edge delamination in curved (0(4)/+/- 45(6))(s) glass-fibre/epoxy beams loaded in bending. Compos Sci Technol 1998; 58: 155–161.

29.

Dugdale

. Yielding of steel sheets containing slits. J Mech Phy Solids 1960; 8(2): 100–104.

30.

Comanho

Dávila

de Moura

. Numerical simulation of mixed-modeprogressive delamination in composite materials. J Compos Mater 2003; 35: 641–652.

31.

Robinson

Galvanetto

Tumino

. Numerical simulation of fatigue-driven delamination using interface elements. Int J Numer Meth Eng 2005; 63: 1824–1848.

32.

Fan

Ben

Cheng

JJR

. Cohesive zone with continuum damage properties for simulation of delamination development in fiber composites and failure of adhesive joints. Eng Fract Mech 2008; 75: 3866–3880.

33.

Motamedi

Mohammadi

. Dynamic analysis of fixed cracks in composites by the extended finite element method. Eng Fract Mech 2010; 77: 3373–3393.

34.

Ashari

Mohammadi

. Modeling delamination in composite laminates using XFEM by new orthotropic enrichment functions. IOP C Ser: Mater Sci Eng 2010; 10: 012240.

35.

Wang

Waisman

. Progressive delamination analysis of composite materials using XFEM and a discrete damage zone model. Comput Mech 2014; 55: 1–26.

36.

Cho

Kim

. Higher-order zig-zag theory for laminated composites with multiple delaminations. J Appl Mech 2001; 68: 869–877.

37.

Cho

. Efficient higher order composite plate theory for general lamination configurations. AIAA J 1993; 31: 1299–1306.

38.

Heidari-Rarani

Sayedain

. Finite element modeling strategies for 2D and 3D delamination propagation in composite DCB specimens using VCCT, CZM and XFEM approaches. Theo Appl Frac Mech 2019; 103: 102246.

39.

Heidari-Rarani

Shokrieh

Camanho

. Finite element modeling of mode I delamination growth in laminated DCB specimens with R-curve effects. Compos Part B: Eng 2013; 45: 897–903.