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.

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{array}{c}\multicolumn{1}{c}{u_1(x_1,x_2,z)=u_1^0(x_1,x_2)+{\psi }_1(x_1,x_2)z}\\ +{\xi }_1(x_1,x_2)z^2+{\varphi }_1(x_1,x_2)z^3\hfill \\ +\sum _{k=1}^{N-1}{S}_1^k(x_1,x_2)(z-z_k)H(z-z_k)\hfill \\ +\sum _{k=1}^{N-1}{\overline{u}}_1^k(x_1,x_2)H(z-z_k)\hfill \end{array}$$

(1)

$$\begin{array}{c}\multicolumn{1}{c}{u_2(x_1,x_2,z)=u_2^0(x_1,x_2)+{\psi }_2(x_1,x_2)z}\\ +{\xi }_2(x_1,x_2)z^2+{\varphi }_2(x_1,x_2)z^3\hfill \\ +\sum _{k=1}^{N-1}{S}_2^k(x_1,x_2)(z-z_k)H(z-z_k)\hfill \\ +\sum _{k=1}^{N-1}{\overline{u}}_2^k(x_1,x_2)H(z-z_k)\hfill \end{array}$$

(2)

$$u_3(x_1,x_2,z)=w(x_1,x_2)+\sum _{k=1}^{N-1}{\overline{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, ${\psi }_i(x_1,x_2)$ indicates the rotation about $x_i$ axis, ${\xi }_i$ and ${\varphi }_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)\mathrm{th}$ layer compared with that of the kth layer at the kth interface while ${\overline{u}}_i^k(x_1,x_2)$ and ${\overline{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, ${\sigma }_{i3}{|}_{z=0}=0$ and ${\sigma }_{i3}{|}_{z=h}=0$ , the transverse strain can be deduced as follow

$$\begin{array}{c}\multicolumn{1}{c}{{\gamma }_{i3}=\frac{\partial u_i}{\partial z}+\frac{\partial u_3}{\partial x_i}={\psi }_i+w_{,i}+2{\xi }_{i}z+3{\varphi }_i{z}^2}\\ +\sum _{k=1}^{N-1}(S_i^k+{\overline{w}}_{,i}^k)H(z-z_k)\hfill \end{array}$$

(4)

On the two surfaces

$${\gamma }_{i3}{|}_{z=0}={\psi }_i+w_{,i}=0$$

(5)

$$\begin{array}{c}\multicolumn{1}{c}{{\gamma }_{i3}{|}_{z=h}={\psi }_i+w_{,i}+2{\xi }_{i}h+3{\varphi }_i{h}^2}\\ +\sum _{k=1}^{N-1}(S_i^k+{\overline{w}}_{,i}^k)=0\hfill \end{array}$$

(6)

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

$${\tau }_{i3}{|}_{z=z_k^{-}}={\tau }_{i3}{|}_{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_{\mathrm{ij}}^k{\varphi }_j-{\overline{w}}_{,i}^k$$

(8)

where $a_{\mathrm{ij}}^k$ is only related to material properties. The derivation of $a_{\mathrm{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{array}{c}\multicolumn{1}{c}{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}{\overline{U}}_i^{k}H(z-z_k)\hfill \end{array}$$

(9)

$$u_3=w+\sum _{k=1}^{N-1}{\overline{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

$$\left\{\begin{array}{c}{\sigma }_{11}^{\left(k\right)}\hfill \\ \multicolumn{1}{c}{{\sigma }_{22}^{\left(k\right)}}\\ \multicolumn{1}{c}{{\sigma }_{23}^{\left(k\right)}}\\ \multicolumn{1}{c}{{\sigma }_{13}^{\left(k\right)}}\\ \multicolumn{1}{c}{{\sigma }_{12}^{\left(k\right)}}\end{array}\right\}=\left[\begin{array}{ccccc}{\overline{Q}}_{1111}^{\left(k\right)}\hfill & {\overline{Q}}_{1222}^{\left(k\right)}\hfill & 0\hfill & 0\hfill & {\overline{Q}}_{1612}^{\left(k\right)}\hfill \\ \multicolumn{1}{c}{{\overline{Q}}_{2211}^{\left(k\right)}}& {\overline{Q}}_{2222}^{\left(k\right)}\hfill & 0\hfill & 0\hfill & {\overline{Q}}_{2212}^{\left(k\right)}\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & {\overline{Q}}_{2323}^{\left(k\right)}\hfill & {\overline{Q}}_{2313}^{\left(k\right)}\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & {\overline{Q}}_{1323}^{\left(k\right)}\hfill & {\overline{Q}}_{1313}^{\left(k\right)}\hfill & 0\hfill \\ \multicolumn{1}{c}{{\overline{Q}}_{1211}^{\left(k\right)}}& {\overline{Q}}_{1222}^{\left(k\right)}\hfill & 0\hfill & 0\hfill & {\overline{Q}}_{1212}^{\left(k\right)}\hfill \end{array}\right]\left\{\begin{array}{c}{\epsilon }_{11}^{\left(k\right)}\hfill \\ \multicolumn{1}{c}{{\epsilon }_{22}^{\left(k\right)}}\\ \multicolumn{1}{c}{{\epsilon }_{23}^{\left(k\right)}}\\ \multicolumn{1}{c}{{\epsilon }_{13}^{\left(k\right)}}\\ \multicolumn{1}{c}{{\epsilon }_{12}^{\left(k\right)}}\end{array}\right\}$$

(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{array}{c}{\sigma }_{\mathrm{ij}}^{\left(k\right)}={\overline{Q}}_{\mathrm{ijmn}}^{\left(k\right)}\hfill \\ [1,z,z^2,z^3,(z-z_p)H(z-z_p),H(z-z_p)]\left\{\begin{array}{c}{\epsilon }_{\mathrm{mn}}^{(0)}\hfill \\ \multicolumn{1}{c}{{\epsilon }_{\mathrm{mn}}^{(1)}}\\ \multicolumn{1}{c}{{\epsilon }_{\mathrm{mn}}^{(2)}}\\ \multicolumn{1}{c}{{\epsilon }_{\mathrm{mn}}^{(3)}}\\ \multicolumn{1}{c}{{\epsilon }_{\mathrm{mn}}^p}\\ \multicolumn{1}{c}{{\overline{\epsilon }}_{\mathrm{mn}}^p}\end{array}\right\}\end{array}$$

(12)

$${\sigma }_{i3}^{\left(k\right)}={\overline{Q}}_{i3m3}^{\left(k\right)}[z,z^2,H(z-z_p)]\left\{\begin{array}{c}{\gamma }_{m3}^{(1)}\hfill \\ \multicolumn{1}{c}{{\gamma }_{m3}^{(2)}}\\ \multicolumn{1}{c}{{\gamma }_{m3}^p}\end{array}\right\}$$

(13)

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

$$\begin{array}{c}\multicolumn{1}{c}{\left\{\begin{array}{c}{N}_{\mathrm{ij}}\hfill \\ \multicolumn{1}{c}{M_{\mathrm{ij}}}\\ \multicolumn{1}{c}{R_{\mathrm{ij}}^{\left(2\right)}}\\ \multicolumn{1}{c}{\begin{array}{c}{R}_{\mathrm{ij}}^{\left(3\right)}\hfill \\ \multicolumn{1}{c}{{\overline{N}}_{\mathrm{ij}}^l}\\ \multicolumn{1}{c}{{\overline{M}}_{\mathrm{ij}}^l}\end{array}}\end{array}\right\}=\left[\begin{array}{cccccc}{A}_{\mathrm{ijmn}}^{\left(0\right)}\hfill & A_{\mathrm{ijmn}}^{\left(1\right)}\hfill & A_{\mathrm{ijmn}}^{\left(2\right)}\hfill & A_{\mathrm{ijmn}}^{\left(3\right)}\hfill & B_{\mathrm{ijmn}}^{k\left(0\right)}\hfill & E_{\mathrm{ijmn}}^{k\left(0\right)}\hfill \\ \multicolumn{1}{c}{A_{\mathrm{ijmn}}^{\left(1\right)}}& A_{\mathrm{ijmn}}^{\left(2\right)}\hfill & A_{\mathrm{ijmn}}^{\left(3\right)}\hfill & A_{\mathrm{ijmn}}^{\left(4\right)}\hfill & B_{\mathrm{ijmn}}^{k\left(1\right)}\hfill & E_{\mathrm{ijmn}}^{k\left(1\right)}\hfill \\ \multicolumn{1}{c}{A_{\mathrm{ijmn}}^{\left(2\right)}}& A_{\mathrm{ijmn}}^{\left(3\right)}\hfill & A_{\mathrm{ijmn}}^{\left(4\right)}\hfill & A_{\mathrm{ijmn}}^{\left(5\right)}\hfill & B_{\mathrm{ijmn}}^{k\left(2\right)}\hfill & E_{\mathrm{ijmn}}^{k\left(2\right)}\hfill \\ \multicolumn{1}{c}{A_{\mathrm{ijmn}}^{\left(3\right)}}& A_{\mathrm{ijmn}}^{\left(4\right)}\hfill & A_{\mathrm{ijmn}}^{\left(5\right)}\hfill & A_{\mathrm{ijmn}}^{\left(6\right)}\hfill & B_{\mathrm{ijmn}}^{k\left(3\right)}\hfill & E_{\mathrm{ijmn}}^{k\left(3\right)}\hfill \\ \multicolumn{1}{c}{B_{\mathrm{ijmn}}^{l\left(0\right)}}& B_{\mathrm{ijmn}}^{l\left(1\right)}\hfill & B_{\mathrm{ijmn}}^{l\left(2\right)}\hfill & B_{\mathrm{ijmn}}^{l\left(3\right)}\hfill & D_{\mathrm{ijmn}}^{\mathrm{lk}}\hfill & F_{\mathrm{ijmn}}^{\mathrm{lk}}\hfill \\ \multicolumn{1}{c}{E_{\mathrm{ijmn}}^{l\left(0\right)}}& E_{\mathrm{ijmn}}^{l\left(1\right)}\hfill & E_{\mathrm{ijmn}}^{l\left(2\right)}\hfill & E_{\mathrm{ijmn}}^{l\left(3\right)}\hfill & F_{\mathrm{ijmn}}^{\mathrm{lk}}\hfill & E_{\mathrm{ijmn}}^{\mathrm{lk}}\hfill \end{array}\right]}\\ \left\{\begin{array}{c}{\epsilon }_{\mathrm{mn}}^{(0)}\hfill \\ \multicolumn{1}{c}{{\epsilon }_{\mathrm{mn}}^{(1)}}\\ \multicolumn{1}{c}{{\epsilon }_{\mathrm{mn}}^{(2)}}\\ \multicolumn{1}{c}{{\epsilon }_{\mathrm{mn}}^{(3)}}\\ \multicolumn{1}{c}{{\epsilon }_{\mathrm{mn}}^k}\\ \multicolumn{1}{c}{{\overline{\epsilon }}_{\mathrm{mn}}^k}\end{array}\right\}\hfill \end{array}$$

(14)

$$\left\{\begin{array}{c}{V}_i^{\left(1\right)}\hfill \\ \multicolumn{1}{c}{V_i^{\left(2\right)}}\\ \multicolumn{1}{c}{Q_i^l}\end{array}\right\}=\left[\begin{array}{ccc}{A}_{i3m3}^{\left(2\right)}\hfill & A_{i3m3}^{\left(3\right)}\hfill & E_{i3m3}^{k\left(1\right)}\hfill \\ \multicolumn{1}{c}{A_{i3m3}^{\left(3\right)}}& A_{i3m3}^{\left(4\right)}\hfill & E_{i3m3}^{k\left(2\right)}\hfill \\ \multicolumn{1}{c}{E_{i3m3}^{l\left(1\right)}}& E_{i3m3}^{l\left(2\right)}\hfill & E_{i3m3}^{\mathrm{lk}}\hfill \end{array}\right]\left\{\begin{array}{c}{\gamma }_{m3}^{\left(1\right)}\hfill \\ \multicolumn{1}{c}{{\gamma }_{m3}^{\left(2\right)}}\\ \multicolumn{1}{c}{{\gamma }_{m3}^k}\end{array}\right\}$$

(15)

where $N_{\mathrm{ij}}$ , $M_{\mathrm{ij}}$ , $R_{\mathrm{ij}}^{\left(2\right)}$ , $R_{\mathrm{ij}}^{\left(3\right)}$ , ${\overline{N}}_{\mathrm{ij}}^l$ , ${\overline{M}}_{\mathrm{ij}}^l$ , $V_i^{\left(1\right)}$ , $V_i^{\left(2\right)}$ , 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

$$\left\{u_i\right\}={\{u_1^0,u_2^0,w,w_{,1},w_{,2},{\varphi }_1,{\varphi }_2,{\overline{u}}_1^1,\dots ,{\overline{u}}_1^{N-1},{\overline{u}}_2^1,\dots ,{\overline{u}}_2^{N-1},{\overline{w}}^k,\dots ,{\overline{w}}^{N-1},{\overline{w}}_{,1}^1,\dots ,{\overline{w}}_{,1}^{N-1},{\overline{w}}_{,2}^1,\dots ,{\overline{w}}_{,2}^{N-1}\}}^T$$

(16)

where $i$ represents the $i\mathrm{th}$ node of the element.

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

$$\begin{array}{c}[u_i^0,{\varphi }_i,{\overline{u}}_i^k]=\sum _{m=1}^n{N}_m[{\left(u_i^0\right)}_m,{\left({\varphi }_i\right)}_m,{\left({\overline{u}}_i^k\right)}_m]\hfill \\ w=\sum _{m=1}^n[P_m{\left(w\right)}_m+H_{\mathrm{xm}}{\left(w_{,x}\right)}_m+H_{\mathrm{ym}}{\left(w_{,y}\right)}_m]\hfill \\ {\overline{w}}^i=\sum _{m=1}^n[P_m{\left({\overline{w}}^i\right)}_m+H_{\mathrm{xm}}{\left({\overline{w}}_{,x}^i\right)}_m+H_{\mathrm{ym}}{\left({\overline{w}}_{,y}^i\right)}_m]\hfill \end{array}$$

(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_{\mathrm{ym}}$ are Hermite interpolation functions. Strain vectors are calculated from nodal displacements of element

$$\begin{array}{c}\multicolumn{1}{c}{\left\{{\epsilon }_{\mathrm{ij}}\right\}={\left[B\right]}_b\left\{u_n\right\}=\sum _{m=1}^n\left[B_{\mathrm{mb},\mathrm{ij}}\right]\left\{u_m\right\}}\\ =[B_{1b,\mathrm{ij}},B_{2b,\mathrm{ij}},\dots ,B_{\mathrm{nb},\mathrm{ij}}]\left\{\begin{array}{c}\left\{u_1\right\}\hfill \\ \multicolumn{1}{c}{\left\{u_2\right\}}\\ \multicolumn{1}{c}{⋮}\\ \multicolumn{1}{c}{\left\{u_n\right\}}\end{array}\right\}\hfill \end{array}$$

(18)

$$\begin{array}{c}\multicolumn{1}{c}{\left\{{\gamma }_{i3}\right\}={\left[B\right]}_s\left\{u_n\right\}=\sum _{m=1}^n\left[B_{\mathrm{ms},\mathrm{ij}}\right]\left\{u_m\right\}}\\ =[B_{1s,\mathrm{ij}},B_{2s,\mathrm{ij}},\dots ,B_{\mathrm{ns},\mathrm{ij}}]\left\{\begin{array}{c}\left\{u_1\right\}\hfill \\ \multicolumn{1}{c}{\left\{u_2\right\}}\\ \multicolumn{1}{c}{⋮}\\ \multicolumn{1}{c}{\left\{u_n\right\}}\end{array}\right\}\hfill \end{array}$$

(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_{\mathrm{ft}}={\left(\frac{{\sigma }_{11}}{X_T}\right)}^2+\frac{{\sigma }_{12}^2+{\sigma }_{13}^2}{S_{12}^2}\ge 1({\sigma }_{11}\ge 0)$$

(20)

$$f_{\mathrm{fc}}={\left(\frac{{\sigma }_{11}}{X_C}\right)}^2\ge 1({\sigma }_{11}<0)$$

(21)

$$\begin{array}{c}\multicolumn{1}{c}{f_{\mathrm{mt}}={\left(\frac{{\sigma }_{22}+{\sigma }_{33}}{Y_T}\right)}^2+\frac{{\sigma }_{23}^2-{\sigma }_{22}{\sigma }_{33}}{S_{23}^2}}\\ +\frac{{\sigma }_{12}^2+{\sigma }_{13}^2}{S_{12}^2}\ge 1({\sigma }_{22}+{\sigma }_{33}\ge 0)\hfill \end{array}$$

(22)

$$\begin{array}{c}\multicolumn{1}{c}{f_{\mathrm{mc}}=({\left(\frac{Y_C}{2S_{23}}\right)}^2-1)\left(\frac{{\sigma }_{22}+{\sigma }_{33}}{Y_C}\right)}\\ +{\left(\frac{{\sigma }_{22}+{\sigma }_{33}}{2S_{23}}\right)}^2+\frac{{\sigma }_{23}^2-{\sigma }_{22}{\sigma }_{33}}{S_{23}^2}\hfill \\ +\frac{{\sigma }_{12}^2+{\sigma }_{13}^2}{S_{12}^2}\ge 1({\sigma }_{22}+{\sigma }_{33}<0)\hfill \end{array}$$

(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:

$$\alpha =\frac{1}{2}{\arctan }^{-1}\frac{2{\tau }_{23}}{{\sigma }_{22}-{\sigma }_{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}$
$\mathrm{ft}$	0.01	–	–	0.01	–	0.01	0.01	–	0.01
$\mathrm{fc}$	0.2	–	–	0.01	0.01	0.01	0.01	–	0.01
$\mathrm{mt}$	–	0.01	–	0.01	–	0.01	–	–	–
$\mathrm{mc}$	–	0.2	–	0.2	0.2	–	0.2	–	–

Inter-laminate damage modes

The delamination failure criteria is applied as equation (25)

$$f_{\mathrm{df}}={\left(\frac{{\sigma }_{33}}{Z_T^{\mathrm{int}}}\right)}^2+{\left(\frac{{\sigma }_{13}}{S_{13}^{\mathrm{int}}}\right)}^2+{\left(\frac{{\sigma }_{23}}{S_{23}^{\mathrm{int}}}\right)}^2\ge 1({\sigma }_{22}+{\sigma }_{33}\ge 0)$$

(25)

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

$${\sigma }_{13}=G_{13}{\gamma }_{13},{\sigma }_{13}=G_{13}{\gamma }_{13}$$

(26)

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

$${\sigma }_{33}=E_2\frac{{\overline{w}}^k}{h_{\mathrm{ply}}}$$

(28)

where $S_e$ is the area of the plate element, and $h_{\mathrm{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 ${\overline{w}}^k$ of all nodes are released in order to calculate ${\sigma }_{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 ${\overline{u}}_i^k$ , ${\overline{w}}_i^k$ , and ${\overline{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 ${\overline{u}}_i^k={\overline{w}}_i^k={\overline{w}}_{,i}^k=0$ . If the delamination occurs between layers, the constrains on ${\overline{u}}_i^k$ , ${\overline{w}}_i^k$ , and ${\overline{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\mathrm{mm}\times 25\mathrm{mm}\times 4.2\mathrm{mm}$ , and a pre-crack with length of $35\mathrm{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\mathrm{mm}\times 1\mathrm{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$ $(\mathrm{GPa})$	$E_2$ $(\mathrm{GPa})$	$G_{12}$ $(\mathrm{GPa})$	$v_{12}$	$Z_T^{\mathrm{int}}$ $(\mathrm{MPa})$	$S_{13}^{\mathrm{int}}$ $(\mathrm{MPa})$	$S_{23}^{\mathrm{int}}$ $(\mathrm{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$ $(\mathrm{GPa})$	$E_2$ $(\mathrm{GPa})$	$E_3$ $(\mathrm{GPa})$	$G_{12}$ $(\mathrm{GPa})$	$G_{13}$ $(\mathrm{GPa})$	$G_{23}$ $(\mathrm{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$ $(\mathrm{MPa})$	$X_C$ $(\mathrm{MPa})$	$Y_T$ $(\mathrm{MPa})$	$Y_C$ $(\mathrm{MPa})$	$Z_T$ $(\mathrm{MPa})$	$S_{12}$ $(\mathrm{MPa})$	$S_{13}$ $(\mathrm{MPa})$	$S_{23}$ $(\mathrm{MPa})$	$G_{\mathrm{IC}}$ $(\mathrm{J}/{\mathrm{m}}^2)$	$G_{\mathrm{IIC}}$ $(\mathrm{J}/{\mathrm{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

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°).