3.3.1 Intra-laminar failure

The criteria used in Abaqus/Explicit are strength-based failure methods, which define activation functions for each of the failure models involved, including fibre, matrix, and loading directions. No dissipation of energy is caused by these criteria. Instead, the functions account for the activation of the yielding functions (>1) describing the failure ‘density’ based on the effective (undamaged) stress field. These predictions are set to reproduce the impact test up to the onset of damage at a pre-existing crack embedded within the seventh specimen interface.

The Hashin, Puck, and Cuntze criteria are used for the analysis of the intra-laminar failure [Reference Hinton, Soden and Kaddour32]. These criteria allow for predicting the failure mechanisms caused by two main modes: Fibre Failure (FF) and Inter-Fibre Failure (IFF). They have been continuously developed since 1980 when Hashin established the first approach [Reference Puck and Schürmann40]. The criteria by Puck [Reference Puck and Schürmann40] and Cuntze [Reference Cuntze and Freund41] are more recent, in which efforts have been made to improve the IFF mode. IFF is important for the current study as it typically appears as part of bending events [Reference Cuntze and Freund41], caused by the impact. The Hashin criterion is based on the modifications of the Mohr/Coulomb theory through the quadratic functions. Equations (1) and (2) show the formulation for the two inter-fibres modes in tension (referred to as ${{IF}}{{{F}}_{{ + }}}$ ) and compression (referred to as ${{IF}}{{{F}}_{{ - }}}$ ), respectively.

(1) \begin{align}IF{F_ + } = \frac{{{{\left( {{\sigma _{22}} + {\sigma _{33}}} \right)}^2}}}{{Y_t^2}} + \frac{{\left( {\sigma _{23}^2 - {\sigma _{22}}{\sigma _{33}}} \right)}}{{S_{23}^2}} + \frac{{\left( {\sigma _{12}^2 + \sigma _{13}^2} \right)}}{{S_{12}^2}}{\rm{\;\;\;\;\;}}for\left( {{\sigma _{22}} + {\sigma _{33}}} \right) \geqslant 0\end{align}

(2) \begin{align}IF{F_ - } & = \frac{1}{{{Y_c}}}\left[ {\left( {\frac{{{Y_c}}}{{2{S_{23}}}}} \right) - 1} \right]\left( {{\sigma _{22}} + {\sigma _{33}}} \right) + \frac{{{{\left( {{\sigma _{22}} + {\sigma _{33}}} \right)}^2}}}{{4S_{23}^2}} + \frac{{\left( {\sigma _{23}^2 - {\sigma _{22}}{\sigma _{33}}} \right)}}{{S_{23}^2}} \nonumber \\[5pt] & \quad + \frac{{\left( {\sigma _{12}^2+ \sigma _{13}^2} \right)}}{{S_{12}^2}}{\rm{\;\;\;\;\;}}for\left( {{\sigma _{22}} + {\sigma _{33}}} \right) \lt 0\end{align}

The Puck criterion improves the IFF mode by calculating it as a function of, the so-called, action plane stresses ( ${\sigma _n},\;{\tau _{nt}},\;{\tau _{nl}}$ ). The common action plane defines the angle of the fracture, ${\rm{\theta }}$ , where the plane stresses act. Unlike Hashin criterion, Puck is not defined with a fixed angle. The fracture angle is calculated by means of an iterative process related to the IFF functions between −90 $^{\circ}$ and 90 $^{\circ}$ that maximises the IFF mode [Reference Puck and Schürmann40]. The IFF mode, as a function of the fracture angle in tensile loading is named ${{IF}}{{{F}}_{{ + }}}\left( {{\theta }} \right)$ , and in compressive loading ${{IF}}{{{F}}_{{ - }}}\left( {\rm{\theta }} \right)$ for the Puck criterion, derivations for which can be found in Ref. (Reference Puck and Schürmann40). The action plane stresses used in the Puck criterion are calculated using Equations (3), (4) and (5) [Reference Puck and Schürmann40].

(3) \begin{align}{\sigma _n}\left( \theta \right) = {\sigma _{22}}{\cos ^2}\theta + {\sigma _{33}}{\rm{si}}{{\rm{n}}^2}\theta + 2{\sigma _{23}}{\rm{sin}}\theta {\rm{cos}}\theta \end{align}

(4) \begin{align}{\tau _{nt}}\left( \theta \right) = \left( {{\sigma _{33}} - {\sigma _{22}}} \right){\rm{sin}}\theta {\rm{cos}}\theta + {\sigma _{23}}\left( {{\rm{co}}{{\rm{s}}^2}\theta - {\rm{si}}{{\rm{n}}^2}\theta } \right)\end{align}

(5) \begin{align}{\tau _{nl}}\left( \theta \right) = {\sigma _{31}}\sin \theta + {\sigma _{21}}{\rm{cos}}\theta \end{align}

The Cuntze criterion is based on the definition of five different failure mechanisms [Reference Cuntze and Freund41, Reference Rezasefat, Gonzalez-jimenez, Giglio and Manes42]. Two of them are related to the fibres, one for tension failure (FF1) and one for compression failure (FF2), with the remaining three associated with the inter-fibre failure, corresponding to the matrix rupture. These include tensile matrix failure (IFF1), compressive matrix failure (IFF2), and matrix shear failure (IFF3). The Cuntze criterion considers the interaction of these modes by using a global failure criterion ( ${{Ef{f}}}{^m}$ ), where each failure mode represents a portion of the load-carrying capacity of the material. These capacities quantify the ‘risk’ of fracture in their respective directions, as well as the risk of the global failure, which is calculated using Equation (6) [Reference Cuntze and Freund41].

(6) \begin{align}Ef{f^m} = {\left( {Ef{f^{\parallel \sigma }}} \right)^m} + {\left( {Ef{f^{\parallel \tau }}} \right)^m} + {\left( {Ef{f^{ \bot \sigma }}} \right)^m} + {\left( {Ef{f^{ \bot \tau }}} \right)^m} + {\left( {Ef{f^{ \bot \parallel }}} \right)^m}\end{align}

The variable ${\rm{m}}$ here is the interaction exponent obtained by curve fitting (for CFRP, the recommended range is 2.5< ${\rm{m}}$ <3.5). The interaction between FF and IFF modes is defined using stress effort, which sums up the contribution of each model. Therefore, the Cuntze’s criterion was defined for the IFF depending on the loading direction as ${{IF}}{{{F}}_{{ + }}}$ and ${{IF}}{{{F}}_{{ - }}}$ , and with a shear stress influence as ${{IF}}{{{F}}_{{{shear}}}}$ . A more detailed description and derivation of Cuntze’s functions can be found in Ref. (Reference Cuntze and Freund41).

In the current study, the same failure modes in the fibre direction (tension and compression) were assumed for the Puck and Cuntze criteria, as shown in Equations (7) and (8).

(7) \begin{align}FF1 = \frac{{{\sigma _{11}}}}{{{X_t}}}\end{align}

(8) \begin{align}FF2 = \frac{{ - {\sigma _{11}}}}{{{X_c}}}\end{align}

The Hashin criterion uses similar failure mode in compression as denoted in Equation (8) above. The tensile fibre mode is different, and is shown in Equation (9):

(9) \begin{align}FF1 = {\left( {\frac{{{\sigma _{11}}}}{{{X_t}}}} \right)^2} + \frac{{{\sigma _{12}}^2 + {\sigma _{13}}^2}}{{{S_{12}}^2}}\end{align}

3.3.2 Inter-laminar failure

VCCT is an efficient method for delamination analysis. The theoretical background of the model is based on the work done by Rybicki and Kanninen [Reference Rybicki and Kanninen43]. The method has been implemented in several software codes, including Abaqus/Explicit. VCCT evaluates the energy release rate ( $G$ ) for each fracture mode at the delamination edge (crack-tip) by using nodal displacements and forces. The criticality of the delamination initiation is defined with the mixed-mode fracture criterion. The mixed-mode fracture criterion is employed in the current work as Equation (10):

(10) \begin{align}f = {\left( {{G_{\rm{I}}}/{G_{{\rm{Ic}}}}} \right)^{\rm{\alpha }}} + {\left( {{G_{{\rm{II}}}}/{G_{{\rm{IIc}}}}} \right)^{\rm{\beta }}}\end{align}

Based on previous investigations with the same material [Reference Hintikka, Wallin and Saarela44], the fracture toughness was defined for Mode I, ${{{G}}_{{{Ic}}}}$ , as 72J/ ${{\rm{m}}^{\rm{2}}}$ and Mode II, ${{{G}}_{{{IIc}}}}{{,}}$ as 779J/ ${{\rm{m}}^{\rm{2}}}$ . The power values for Mode I ( ${\rm{\alpha }}$ ) and II ( ${\rm{\beta }}$ ) were defined as 0.75 and 0.69, respectively. A high value for G _IIIC (10,000J/m²) is used in analysis to remove the effects of fracture Mode III. VCCT is modelled between the bottom two plies, seventh and eighth. A pre-crack is defined for the interface plane.

3.4.1 Intra-laminar failure

In LS-Dyna, MAT54 material model is used to model the unidirectional laminate. The stress-strain behaviour of the material in the elastic region is assessed by Equations (11)–(13) [Reference Wade, Feraboli, Osborne and Rassaian35]:

(11) \begin{align}{\varepsilon _1} = \frac{1}{{{E_1}}}\left( {{\sigma _1} - {_{12}}{\sigma _2}} \right)\end{align}

(12) \begin{align}{\varepsilon _2} = \frac{1}{{{E_2}}}\left( {{\sigma _2} - {v_{21}}{\sigma _1}} \right)\end{align}

(13) \begin{align}2{\varepsilon _{12}} = \frac{1}{{{G_{12}}}}{\tau _{12}} + c{\tau _{12}}^3\end{align}

where direction-1 indicates the fibre axial direction, direction-2 indicates the matrix transverse direction, and direction-12 indicates the shear direction. Here, $c$ is an input parameter accounting for the nonlinear shear stress term, which must be calibrated whenever shear is present.

The Chang-Chang composite damage model [Reference Chang and Kuo-yen34] is used within MAT54 beyond the elastic region. It is an updated form of Hashin comprehensive composite damage model [Reference Hashin45]. The model assumes that unidirectional CFRP behaves transversely isotropic in the fibre direction. It includes a mixed mode term, known as the shear stress weighing factor, ${\rm{\beta }}$ , which allows interaction between shear and normal failure modes. Setting ${\rm{\beta }}$ to 0 initiates the Hashin failure criterion, whereas setting it to 1 executes maximum stress failure criteria [Reference Wade, Feraboli, Osborne and Rassaian35]. This damage model analytically computes the stresses and damage progression of the fibre and matrix separately to determine the condition of an individual ply. The Chang-Chang failure mode is described by Equations (14)–(17) where the tensile fibre mode is as follows:

(14) \begin{align}{e_f}^2 = {\left( {\frac{{{\sigma _{11}}}}{{{X_t}}}} \right)^2} + \beta {\left( {\frac{{{\sigma _{12}}}}{{{S_c}}}} \right)^2} - 1{\rm{\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;}}\left\{ {\begin{array}{*{20}{c}}{{e_f}^2\;\; \geqslant \;\;0\; \to \;failed\;}\\{{e_f}^2\; \lt \;0\; \to \;elastic}\end{array}} \right.\end{align}

The compressive fibre mode is defined as:

(15) \begin{align}{\rm{\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;}}{e_c}^2 = {\left( {\frac{{{\sigma _{11}}}}{{{X_c}}}} \right)^2} - 1{\rm{\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;}}\left\{ {\begin{array}{*{20}{c}}{{e_c}^2\;\; \geqslant \;\;0\; \to \;failed\;}\\{{e_c}^2\; \lt \;0\; \to \;elastic}\end{array}} \right.\end{align}

where ${{{X}}_{{t}}}$ and ${{{X}}_{{c}}}$ are the fibre longitudinal tensile and compressive strengths, respectively, and ${{{e}}_{{f}}}$ as well as ${{{e}}_{{c}}}$ are the fibre failure parameter. ${{\rm{\sigma }}_{{\rm{11}}}}$ represents the longitudinal stress and ${{\rm{\sigma }}_{{\rm{12}}}}$ the shear stress of each layer. S _c denotes the matrix shear strength [Reference Hallquist33, Reference Brandon46]. The tensile matrix mode is defined as:

(16) \begin{align}{\rm{\;\;\;\;\;\;\;\;\;\;}}{e_m}^2 = {\left( {\frac{{{\sigma _{22}}}}{{{Y_t}}}} \right)^2} + {\left( {\frac{{{\sigma _{12}}}}{{{S_c}}}} \right)^2} - 1{\rm{\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;}}\left\{ {\begin{array}{*{20}{c}}{{e_m}^2\;\; \geqslant \;\;0\; \to \;failed\;}\\{{e_m}^2\; \lt \;0\; \to \;elastic}\end{array}} \right.\end{align}

And the compressive matrix mode is as follows:

(17) \begin{align}{e_d}^2 = {\left( {\frac{{{\sigma _{22}}}}{{2{S_c}}}} \right)^2} + \left[ {{{\left( {\frac{{{Y_c}}}{{2{S_c}}}} \right)}^2} - 1} \right]\frac{{{\sigma _{22}}}}{{{Y_c}}} + {\left( {\frac{{{\sigma _{12}}}}{{{S_c}}}} \right)^2} - 1{\rm{\;\;\;}}\left\{ {\begin{array}{*{20}{c}}{{e_d}^2\;\; \geqslant \;\;0\; \to \;failed\;}\\{{e_d}^2\; \lt \;0\; \to \;elastic}\end{array}} \right.\end{align}

where ${{{Y}}_{{t}}}$ and ${{{Y}}_{{c}}}$ are the transverse tensile and compressive strength of the matrix. ${{{e}}_{{m}}}$ and ${{{e}}_{{d}}}$ are the matrix failure parameter that determines matrix cracking [Reference Hallquist33–Reference Wade, Feraboli, Osborne and Rassaian35].

The Chang-Chang composite damage criteria within MAT54 material model in LS-Dyna assumes a linear elastic ply level response until failure, without pre- or post-peak softening. It contains a set of non-physical input parameters which determines the element failure and can be divided into three categories: erosion, crashfront softening factors and parameters describing the material behaviour after failure initiation [Reference Aleksandr, Montesano and Butcher47]. For solid elements, an element is deleted when its single integration point has met the failure criteria. Elements with shared nodes with deleted element become crashfront elements which can have their strengths reduced by varying the non-physical SOFT input parameters. Failed ply can still carry a significant amount of stress and energy, and only once it reaches failure using one of the non-physical input parameters or an effective failure strain, the stresses are reduced to zero and the element is deleted [Reference Hallquist33]. For this study, the built-in default YCFAC strength reduction factor for compressive fibre strength after matrix failure is used to degrade the fibre strengths of the plies if compressive matrix failure would take place [Reference Hallquist33]. This parameter simulates damage caused to the fibres from failed matrix [Reference Wade, Feraboli, Osborne and Rassaian35]. Additionally, the built-in default SOFT material strength reduction parameter is used to reduce the strengths of the elements ahead of the crashfront [Reference Hallquist33]. It is used to degrade the strengths of the surrounding elements to simulate the damage propagation from the crashfront [Reference Wade, Feraboli, Osborne and Rassaian35]. The failure stress and strain inputs were not defined within the composite damage model to ensure an elastic response of the plies during the impact.

3.4.2 Inter-laminar failure

3.4.2.1 Cohesive zone method

CZM is based on the fracture model originally proposed by Dugdale [Reference Dugdale48] and Barenblatt [Reference Barenblatt49]. An alternative to the VCCT, it is one of the most widely used approaches to investigate interface bonding failure [Reference Siddens and Bayandor22–Reference Bayandor, Thomson and Callus24, Reference Brandon46]. It assumes a cohesive damage zone that develops near the crack tip and uses strain energy release rates, ${{G}}$ , during the formation of fracture to predict delamination [Reference Siddens and Bayandor22]. Hence, CZM is governed by the properties of the material, crack initiation condition and a crack evolution function, where it relates the surface loads in normal and shear direction, ${\rm{\sigma }}$ , to the displacement. Crack initiation takes place when the model reaches maximum interface strength, ${{\rm{\sigma }}^{\rm{0}}}$ , that is the maximum load on stress-displacement relation [Reference Siddens and Bayandor22, Reference Turon, Davila, Camanho and Costa50].

Several linear and non-linear constitutive failure damage formulations exist for defining the cohesive law. The law used in this study is governed by a bilinear stress-displacement relation as well as mixed-mode to represent the interactions of Mode I, Mode II, and Mode III separation. Figure 6 shows the bilinear stress-displacement relation used in this study [Reference Siddens and Bayandor22, Reference Bak, Sarrado, Turon and Costa51]. The area under the stress-displacement relation is equal to the energy release rate, ${{G}}$ , or the fracture toughness ${{{G}}_{{c}}}$ , that is the critical value of ${{G}}$ for delamination growth [Reference Bayandor, Thomson and Callus24, 52]. In LS-Dyna, energy release rates are assigned for both Mode I and Mode II. Mode III is assumed to be nearly identical to Mode II separation. This allows the cohesive model to utilise mixed mode capability which combines the energy release rate. This is used in conjunction with a damage formulation (progressive softening) where mixed-mode displacement for total failure is computed using either Power law or Benzeggagh-Kenane (B-K) law [Reference Benzeggagh and Kenane53]. The interface bond softening happens when the separation and sliding between the plies exceeds the softening strain, ${{\rm{\delta }}^{\rm{0}}}{\rm{.\;}}$ Once the failure strain, ${{\rm{\delta }}^{\rm{F}}}$ , is reached, the interface bond fails and delamination occurs [Reference Siddens and Bayandor22]. The formulation for Power and B-K law implemented in LS-Dyna can be found in Ref. (Reference Hallquist33).

Figure 6. Bilinear traction-displacement relationships [Reference Siddens and Bayandor22, Reference Hallquist33].

The slope of the bilinear stress-displacement before damage initiation, K, is known as the interface stiffness or penalty stiffness as it is referred to in LS-Dyna. Turon et al. [Reference Turon, Davila, Camanho and Costa50] formulated that the effective elastic properties of the composite laminate are not affected by the cohesive surface properties whenever ${{{E}}_{{3}}} \ll {{\;Kt}}$ , refer to Equation (18).

(18) \begin{align}K = \frac{{\alpha {E_3}}}{t}\end{align}

where ${{{E}}_{{3}}}$ is the through-thickness Young’s modulus of the material, which for a transversely isotropic material is ${{{E}}_{{3}}}{{ = \;}}{{{E}}_{{2}}}$ and $t$ is the sub-laminate thickness. ${\rm{\alpha }}$ is a parameter that is recommended to be much larger than 1. Values of $\alpha $ are recommended to be greater than 50 [Reference Turon, Davila, Camanho and Costa50].

The crack propagates when the energy release rates reach their critical values ${{{G}}_{{c}}}$ , that is their fracture toughness [Reference Turon, Davila, Camanho and Costa50]. The distance from the crack tip to the point where the maximum stress (cohesive traction) is achieved is referred to as the length of the cohesive zone, ${{{l}}_{{{cz}}}}$ . Several different models have been proposed to estimate the cohesive zone length [Reference Dugdale48, Reference Barenblatt49, Reference Irwin54–Reference Song, Dávila and Rose59]. All of these proposed models have the form as shown in Equation (19).

(19) \begin{align}{l_{cz}} = ME\frac{{{G_c}}}{{{{\left( {{\sigma ^0}} \right)}^2}}}\end{align}

where ${{E}}$ is the Young’s modulus of the laminate and ${{M}}$ is a parameter that depends on the different cohesive zone length models. This study uses the Rice [Reference Rice56] and Falk et al. [Reference Falk, Needleman and Rice57] recommended parameter defining the cohesive zone length. In terms of LS-Dyna parameters, the cohesive zone lengths can be defined for both Mode I and Mode II as shown in Equation (20), where i represents either Mode I or Mode II, and X represents either T (the normal peak load input value) for Mode I, or S (the shear peak load input value) for Mode II.

(20) \begin{align}{l_{cz,i}} = ME\frac{{{G_{Ic}}}}{{{{\left( X \right)}^2}}}\end{align}

For mixed-mode, the cohesive zone length must satisfy the following conditions shown in Equation (21).

(21) \begin{align}{l_e} \le \frac{{{l_{cz,i}}}}{{{N_e}}}\end{align}

where ${{{N}}_{{e}}}$ is the number of elements in the cohesive zone and ${{{l}}_{{e}}}$ is the mesh size used in the direction of crack propagation [Reference Turon, Davila, Camanho and Costa50, Reference Song, Dávila and Rose59]. If the cohesive zone is not discretised by enough elements, then the distribution of load ahead of the crack is not represented correctly in FEM. Therefore, Turon et al. [Reference Turon, Davila, Camanho and Costa50] suggest using a minimum of three elements in the cohesive zone to predict the propagation of delamination. However, such fine mesh requirements can make structural analysis computationally expensive and therefore not feasible for some studies [Reference Song, Dávila and Rose59]. It can be seen from Equation (21) that the cohesive zone length is inversely proportional to the square of the interface peak loads. The length of the cohesive zone can therefore be theoretically lengthened to ensure that they span enough elements of the given size. Using this methodology, we can obtain the required normal and shear peak loads that would allow us to use an appropriate cohesive zone length without making the mesh overly fine [Reference Turon, Davila, Camanho and Costa50, Reference Song, Dávila and Rose59]:

(22) \begin{align}\left\{ {\begin{array}{*{20}{c}}{{T^0} = \;\sqrt {\dfrac{{9\pi E{G_{Ic}}}}{{32{N_e}{l_e}}}} }\\[11pt] {{S^0} = \;\sqrt {\dfrac{{9\pi E{G_{IIc}}}}{{32{N_e}{l_e}}}} }\end{array}} \right.\end{align}

and the normal and shear peak loads to be used in the model can be selected based on the criteria:

(23) \begin{align}\left\{ {\begin{array}{*{20}{c}}{\bar T = Min\left\{ {{T^0},T} \right\}}\\{\bar S = Min\left\{ {{S^0},S} \right\}}\end{array}} \right.\end{align}

3.4.2.2 Tiebreak contact

An alternative to the cohesive zone model technique is the tiebreak contact algorithm in LS-Dyna that uses a segment-based approach to model cohesive interaction between the plies or laminates. By implementing the Dycoss Discrete Crack Method [Reference Lemmen and Meijer60] within the tiebreak contact, a bilinear stress-displacement relation can be used similar to the cohesive zone method. LS-Dyna defines tiebreak contacts as penalty-based with the added feature of including a failure criterion [Reference Hallquist33]. The tiebreak contact works in a similar way as regular contacts under compressive load. The separation phenomenon between surfaces is defined by a strength-based failure criterion [Reference Xiao, Botkin and Johnson61]. For that, it requires definition of a normal failure stress parameter, NFLS, and a shear failure stress parameter, SFLS [Reference Hallquist33, Reference Brandon46].

Tiebreak contact essentially models a theoretical interlaminar bonding based on the normal and shear peak stress. This allows the delamination to propagate from the initiation zone based on mixed-mode, which combines the crack opening damage (mode I) and in-plane shear damage (mode II). This is used in conjunction with a damage formulation where mixed-mode displacement for total failure is computed using either Power Law or Benzeggagh-Kenane’s [Reference Benzeggagh and Kenane53]. An advantage of using tiebreak contact is that it does not require any additional layer of elements to be modelled or defined at the interface as it uses a segment-based approach, where the faces of the ply elements are bonded. Figure 7 shows a comparison between the cohesive zone model and tiebreak contact.

Figure 7. Comparison of cohesive elements and tiebreak contact [Reference Mohammed, Paiva and Mayer62].

4.2.1 Delamination detection

Ultrasonic detection pointed out at a major change in the pulse echo after reaching the boundary condition (the metal clamp) during testing. The location of the change in the observed pulse echo, indicative of the delaminated region, was marked. A circular delamination pattern was identified based on the recorded data. The radius of the delaminated region was measured with the accuracy of ±5mm, as 12.3 ±5mm. Figure 16 shows the specimen with the location of the pulse echo inflection point highlighted.

Figure 16. Inflection point of observed pulse echo identifying delaminated region using 5MHz probe.

4.2.2 Abaqus/Explicit impact model

For the dynamic impact test analysis in Abaqus/Explicit, the initiation of delamination is modelled at the seventh ply interface, where a pre-crack is defined. Figure 17 shows the von Mises stresses through the thickness when delamination starts to evolve. Figure 18 shows the von Mises stress field across the surface of the seventh ply as well as the VCCT results at the seventh ply interface. Results indicate peak stresses near the centre region as expected, where the impactor comes into contact with the plies. The von Mises stresses are presented since they offer a reasonable overall behaviour of the impacted specimen. However, the in-plane stresses (x, y, and z-directions) from the impact simulations are also shown in Fig. 19. The failure initiation occurs in the fibre-matrix planes, so that the in-plane stresses in x- and y-directions are more critical. Results show that the x-plane stresses are the largest.

Figure 17. von Mises stress through thickness when delamination starts to evolve at seventh ply interface.

Figure 18. von Mises stress and VCCT results when delamination starts to evolve at seventh ply interface.

Figure 19. In-plane stresses from VCCT simulation when delamination starts to evolve at seventh ply interface.

Figure 20 shows the activation functions (IFF mode) of the three failure criteria through the thickness and with the pre-crack at the seventh ply interface. The tensile (membrane) response acts at the lower part and compression at the upper part of the laminate due to bending. Figure 21 highlights the failure status (inter- and intra-laminar) of the lower part of the seventh ply. Results are depicted when delamination starts to evolve and include the IFF mode values at the ply. For better comparison, results are normalised to a fixed range from 0 to 1, where 0 denotes “intact” elements and 1 indicates elements that have reached “failure”.

Figure 20. Stress and activation functions status through cross-section of laminate.

Figure 21. Comparison of failure prediction when delamination starts to evolve at seventh ply interface.

In terms of the stress and failure values, the simulated prediction is highlighted by the 48 elements (refer to Fig. 18) which are used for reporting the performances, as provided in Table 4 through Table 8. The simulation is also conducted without the pre-crack (no VCCT) at the interface between seventh and eighth plies using the tie constraint model, results of which are included in the same tables for comparison. The Hashin criterion, shown in Table 5, presents the highest value for the compressive mode, making it the most conservative with respect to the two other criteria. This is due to the definition of the fracture plane at a ‘pre-fixed’ angle (0 $^{\circ}$ wrt. thickness direction), on which transverse stresses act. The simulation without pre-crack shows a slight increase in the probability of tensile failure, as well as an abrupt decrease in the compressive failure due to the stress field reduction (-66.1%).

Table 4. Averaged values of von Mises stresses of highlighted elements at seventh ply interface

Table 5. Hashin criterion, averaged values of IFF mode of marked elements

In the marked elements (refer to Fig. 18), Puck criterion on average results in the formation of a fracture plane at 55.52 $^{\circ}$ (wrt. thickness direction). The results with a pre-existing crack for the compressive mode are found to be lower than Hashin criterion with pre-crack, whereas the results with pre-crack for the tensile mode are found to be higher than Hashin criterion with one (see Table 6). However, without the pre-crack, the IFF values are found to be higher. This solution represents a clear difference in comparison to Hashin criterion.

Table 6. Puck criterion, averaged values of IFF mode of marked elements with calculated fracture angle $\theta_{frac} = 55.622^{\circ}$

The global response (EFF) of Cuntze criterion offers a higher predicted failure value for the case without a pre-crack (see Table 7). The local response of this criterion, unlike Puck, is defined with ‘fixed’ fracture angles based on the experimental results [Reference Cuntze and Freund41]. The IFF1 value is generated by a fracture plane parallel to the fibres (0 $^{\circ}$ wrt. thickness direction) and it reports a slightly higher value in the simulation without a pre-crack. The IFF3 mode is defined by a fracture plane of 53 $^{\circ}$ [Reference Puck and Schürmann40], which results in a higher value than in the case with a pre-crack. The IFF2 (shear) mode is also represented by a fracture plane parallel to the fibres (0 $^{\circ}$ wrt. thickness direction), resulting in the highest simulation values with a pre-crack (i.e., higher stress field). However, the failure value reduces drastically for the simulation without a pre-crack (i.e., lower stress field). Therefore, a similar response to that of the Puck is obtained: that is an increase in the IFF failure prediction without a pre-crack.

Table 7. Cuntze criterion, averaged values of IFF mode of marked elements

The difference between the results of Hashin criterion versus Puck and Cuntze criteria, for the IFF mode, is essentially due to the definition of the fracture angle. Hashin criterion causes a higher IFF value in the simulation with the pre-crack (i.e., higher stress field). However, when the stress field reduces, as in the simulation without a pre-crack, the IFF value decreases. This behaviour does not occur in the other two criteria. However, the two criteria result in an increase in the IFF mode as a consequence of the rise in the failure throughout the UD-plies, rather than the interface, where no delamination takes place (no pre-crack).

As seen in Table 8, for the common responses, including the fibre failure modes (tension/compression) measured by the three criteria, higher values are obtained at the element level when a pre-crack exists. This is due to the partial lack of load-carrying capacity in the plies in the presence of a pre-existing crack. For simulations without a pre-crack, FF modes resulted in values nearing zero. This does not occur for the global response (EFF) of Cuntze criterion, as indicated before (refer to Table 7). For the case of Hashin, the FF tensile mode, defined in Equation (9), results in the same values for pre-crack simulations as those without pre-crack.

Table 8. Averaged values of FF mode (Hashin, Puck, and Cuntze) in marked elements at seventh ply interface

4.2.3 LS-Dyna impact model

The Chang-Chang composite damage model is applied to the plies, and delamination is modelled between the seventh and eighth plies to determine the damage onset. Element deletion is not used to obtain an elastic response. A layer of cohesive elements is modelled between bottom plies with a thickness of 0.01 mm. Unlike VCCT, CZM does not require the existence of an initial crack in the structure and therefore is able to represent the complete interlaminar bonding layer. Figure 22 shows the delamination initiation (a) and the subsequent delamination as it propagates (b). Cohesive elements reaching failure due to the applied load over time are shown, representing the crack initiation and thereby highlighting delamination onset.

Figure 22. (a) Initial delamination at 0.446 ms (b) delamination propagation.

Figure 23 shows the von Mises stresses through the thickness of the laminate and provides a comparison between the CZM and tiebreak models. Figure 24 depicts the comparison of von Mises stresses on the surface of the seventh ply between the CZM and tiebreak constitutive model when delamination initiates. The von Mises stresses are presented since they offer a reasonable overall behaviour of the impacted specimen.

Figure 23. Comparison of von Mises stresses through thickness of modelled specimen between CZM (cohesive elements) and tiebreak model at the seventh ply at delamination onset.

Figure 24. Comparison of the von Mises stresses between CZM (cohesive elements) and tiebreak model at seventh ply at delamination onset.

The stresses have about the same magnitude and their distribution is almost identical. The cohesive zone model predicts slightly higher stresses in the laminate near the centre region as compared to the tiebreak model, however, the difference is minimal. The maximum stresses are found to be at similar locations at which they reach a value of about 1.27 GPa. Near the centre region of the impact, the stresses in both models are found to be in the range of 350–650 MPa. The in-plane stresses (x, y, and z-directions) from the simulated impacts are also provided in Fig. 25. Given that the focus of the study was to investigate the damage in the matrix and fibre action planes, the in-plane stresses in x- and y-direction are more critical. Results show that the x-plane stresses are the largest, similar to the Abaqus/Explicit predictions (refer to Fig. 19).

Figure 25. Comparison of in-plane stresses between CZM (cohesive elements) and tiebreak model at seventh ply at delamination onset.

Figure 26 depicts the fibre and matrix damage when delamination initiates on the bottom side of the seventh ply for the model with cohesive elements. Damage initiation takes place at the centre region where the impactor strikes the laminate as expected. Fibre is shown to be damaged in tension in the + 45 $^{\circ}$ direction along the fibre orientation. The damage further develops near the centre of the fibre in compression. The seventh ply fibres do not seem to fail from a 5J impact since no element reaches the failure value of 1. The matrix on the other hand is seen to be both damaged as well as failed in tension near the centre region of the specimen where the impact primarily takes place. This is expected since in general matrix is more prone to damage and failure than fibre when the ply is deflected downwards due to the applied impact load. In this case, the matrix did undergo damage and failed in compression. The regions where the matrix failed was close to where the delamination was initiated.

Figure 26. Damage at delamination onset.

Figure 27 shows the damage through the thickness of the laminate when delamination initiates for the model with cohesive elements. The damage for fibres in tension starts at the bottom two plies and then propagates upwards. Fibres in compression, on the other hand, are seen to have failed in some regions in the upper plies as expected. The matrix in tension follows a similar trend in which the failure initiates at the bottom plies and propagates upwards, where matrix in compression follows the trend of the fibres in the sense that it fails in the upper plies and undergoes damage in the bottom plies.

Figure 27. Through thickness damage at delamination onset.

Figure 28 shows the force-time history of the CZM and tiebreak constitutive models in LS-Dyna compared with the experimental data during the initial stages of damage. Delamination initiation occurs at around 0.446ms, which is indicated by the black arrow. Both models agree with the experimental force-time history well. Large drops in applied loading are seen around 0.43–0.45 ms when the initial delamination occurs. Further drops in the load histogram indicate delamination propagation and the formation of the initial crack. The simulation results, particularly those of the tiebreak model, do not capture any significant load drops as distinctively observed in the experiment. This may have been caused by the fact that delamination measures were not incorporated into the interfaces between all plies and that the element failure and deletion was not activated in the models. As a result, the laminate can only experience an elastic response, leading to an increased load-bearing capacity of the structure.

Figure 28. Comparison of force-time histories between CZM and tiebreak model in LS-Dyna during the initial stages of delamination.